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Preface 

The  purpose  of  this  thesis  was  to  demonstrate  the  practicality  of 
conducting  an  analysis  of  the  accessibility  of  Earth-approaching  aster¬ 
oids  on  a  micro-computer.  The  work  for  this  thesis  combined  aspects 
of  astronomy  and  computers,  in  which  my  current  interests  are,  respec¬ 
tively,  asteroids  and  the  Turbo  PASCAL  programming  language.  A  reasonably 
ambitious  computer  program  to  integrate  these  interests  into  a  practical 
tool  for  meaningful  analyses  seemed  a  worthwhile  challenge. 

The  computer  program  calculates  the  change  in  velocity  that  must 
be  imparted  to  a  spacecraft  to  transfer  it  from  an  orbit  around  the  Earth 
to  an  orbit  around  an  asteroid.  Although  the  program  will  work  for  any 
asteroid,  only  Earth-approaching  asteroids  were  analyzed.  Particular 
attention  was  given  to  the  three  Earth-approaching  asteroids  that  were 
discovered  from  January  to  November  1985.  Input  to  the  program  is  a 
data  file  containing  the  classical  orbital  elements  of  the  Earth  and  the 
asteroids.  Output  is  via  plots  and  tables  that  show  the  velocity  changes 
required  and  the  launch  opportunities  that  will  be  available  up  to  the 
year  2020  to  achieve  these  intercept  trajectories. 

As  a  relative  novice  to  orbital  mechanics,  I  am  indebted  to  my 
thesis  advisor  Dr.  William  Wiesel  for  his  very  knowledgeable  guidance 
throughout  this  experience.  I  was  fortunate  to  find  an  advisor  with  a 
previous  academic  interest  in  asteroids  and  I  hope  that  I  may  have  some¬ 
what  rekindled  that  interest. 

I  must  also  acknowledge  an  important  contribution  to  my  thesis  from 
Mr.  Neal  Hulkower,  whose  work  is  often  cited  in  this  thesis.  At  a 


ii 


critical  stage  in  this  work,  a  couple  of  long-distance  telephone  con¬ 
versations  with  Mr.  Hulkower  firmly  established  the  way  to  proceed. 

To  my  wife  and  two  young  children  go  my  sincere  appreciation  for 
the  understanding  and  support  throughout  these  last  18  months  of 
academia.  Now  maybe  the  kids  can  play  with  the  computer  again. 

Philip  w.  Somers 
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Abstract 

The  purpose  of  this  paper  was  to  analyze  the  accessibility  of  Earth- 
approaching  asteroids  using  a  computer  program  that  was  practical  to  run 
on  a  micro-computer.  This  analysis  employs  techniques  that  can  easily 
be  adapted  to  find  optimal  trajectories  for  a  variety  of  orbital  inter¬ 
cept  applications. 

The  mathematical  analysis  was  adapted  from  recently-developed 
algorithms  that  were  designed  to  run  on  main  frame  computers  using  exten¬ 
sive  software  libraries  and  data  resources.  The  computer  program 
developed  for  this  paper  was  designed  to  operate  on  an  IBM  PC  equipped 
with  an  8087  math  co-processor  chip.  Programming  was  done  in  Turbo 
PASCAL  Version  3.0  which  supports  the  8087  mathematical  capabilities. 

The  program  was  designed  to  be  self-contained  except  for  data  files  of 
orbital  elements.  The  program  was  also  designed  to  operate  efficiently 
and  quickly  while  retaining  much  of  the  accuracy  found  on  the  main  frame 
implementations.  Only  nonperturbed  Keplerian  motion  was  modelled.  Every 
effort  was  made  to  ensure  the  program  was  as  flexible  as  possible.  Any 
object  in  the  solar  system  in  heliocentric  orbit  can  be  used  as  either 
the  departure  or  arrival  body.  Orbital  element  data  files  are  included 
for  all  the  planets,  several  periodic  comets,  all  the  recently-discovered 
Earth-approaching  asteroids,  and  many  of  the  main  belt  asteroids.  This 
flexibility  permits  not  only  rendezvous  missions  to  be  calculated,  but 
can  just  as  easily  handle  fly-by  trajectories  and  return-to-Earth 
missions. 
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Extensive  testing  of  the  algorithms  was  conducted  on  all  of  the 
Earth-approaching  asteroids  discovered  since  1982,  and  on  several  other 
Earth-approaching  asteroids  that  have  particularly-accessible  orbits. 
Tabular  and  graphical  data  was  formatted  and  displayed  so  as  to  be  highly 
readable  and  practical  to  the  user.  Extensive  use  was  made  of  the  con¬ 
cept  of  Prime  Rib  curves  to  plot  global  optimum  trajectory  changes  of 
velocity  for  all  combinations  of  true  anomalies  of  the  departure  and 
arrival  planets.  Validation  of  the  program  was  done  by  comparison  of  the 
output  to  published  results  from  calculations  done  on  main  frame  com¬ 
puters,  and  by  comparison  of  ephemerides  with  those  published  in  the 
Astronomical  Almanac.  The  results  demonstrate  the  practicality  of 
quickly  and  accurately  analyzing  the  accessibility  of  Earth-approaching 
asteroids  on  a  micro-computer. 


AN  ANALYSIS  OF  THE  ACCESSIBILITY  OF  EARTH-APPROACHING  ASTEROIDS 


I .  Introduction 

A  rendezvous  mission  to  any  of  the  3500  known  asteroids  is  possible, 
but  the  easiest  and  most  practical  missions  would  be  to  the  Earth- 
approaching  asteroids.  However,  no  interplanetary  probes  have  ever 
travelled  to  an  asteroid.  All  of  the  information  currently  known  about 
these  minor  planets  has  come  from  ground-based,  or  near-Earth-orbit 
sensors.  Probes  such  as  Pioneer  and  Voyager  have  gone  to  most  of  the 
planets.  Spacecraft  have  landed  on  Mars  and  Venus,  and  flyby  missions 
have  photographed  Mercury,  Jupiter  and  Saturn.  However,  asteroids  remain 
totally  unexplored.  Most  asteroids  are  much  closer  to  the  Earth  than 
Jupiter  or  Saturn,  and  many  are  easier  to  reach  than  any  of  the  planets. 

In  fact,  some  projected  trips  to  land  on  asteroids  and  return  to  Earth 
would  require  a  smaller  expenditure  of  energy  than  a  similar  trip  to  the 
Moon.  Of  the  82  known  Earth-approaching  asteroids,  28  have  been  dis¬ 
covered  in  the  1980s,  kindling  a  renewed  interest  in  a  rendezvous  mission. 

Algorithms  have  been  developed  to  determine  the  accessibility  (the 
energy  requirements  and  the  frequency  of  launch  opportunities)  of  aster¬ 
oids.  lo  ensure  accuracy  and  timeliness  of  the  calculations,  these 
algorithms  have  been  implemented  on  large  main  frame  computers  with 
extensive  software  and  data  libraries.  A  practical  implementation  is 
needed  that  can  be  run  on  more  modest  computing  facilities.  With  the 
many  recent  discoveries  of  Earth-approaching  asteroids  (10:321-328; 
21:1-19),  and  the  renewed  interest  in  asteroid  exploration,  a 


microcomputer  implementation  would  greatly  improve  the  utility  of  these 
algorithms.  Accessibility  information  would  then  be  available  to  many 
users,  for  the  asteroids  in  which  they  had  an  interest.  This  information 
would  also  be  widely  available  for  newly-discovered  asteroids  as  soon 
as  the  orbital  elements  were  known. 


II.  Background 


Many  reasons  have  been  offered  to  justify  a  rendezvous  with  an 
asteroid.  Some  scientists  believe  that  the  asteroids  are  made  of  very 
primitive  material  largely  undisturbed  since  the  formation  of  the  solar 
system.  Important  scientific  research  could  be  accomplished  with  even  a 
tiny  sample  of  this  material.  Due  to  the  large  number  of  known  aster¬ 
oids,  there  is  expected  to  be  a  correspondingly-large  variety  to  the  com¬ 
position  of  their  material.  According  to  Stancaai  and  Soldner  (18:1), 
the  principle  goal  of  missions  to  asteroids  would  be  to  determine  chemi¬ 
cal  composition  and  structure  to  provide  data  for  the  study  of  the  early 
history  of  the  solar  system.  Due  to  the  extremely  low  gravitational 
forces,  several  factors  that  cause  chemical  changes  in  the  materials  of 
the  large  planets  are  not  a  factor  on  asteroids.  With  no  atmosphere, 
there  is  no  wind  or  rain  erosion  or  surface  decomposition.  Internal 
gravitational  pressures  are  too  low  to  produce  hot  liquid  cores,  except 
possibly  to  some  degree  in  a  few  of  the  largest  asteroids.  The  processes 
of  heating,  melting  and  differentiation  of  materials  and  elements  that 
have  radically  changed  terrestrial  planets  probably  have  not  occurred 
with  asteroids  (23:434).  Small  size  and  low  gravity  make  the  asteroids 
small  targets  for  impacts  by  crater-producing  objects.  However,  they 
have  no  atmosphere  to  protect  them  from  micrometeorites  and  must  suffer 
constant  bombardment.  The  effects  of  this  bombardment  would  probably 
be  confined  to  a  thin  layer  on  the  surface.  Beneath  this  thin  layer, 
the  material  is  expected  to  be  more  volatile  than  that  of  the  inner 
planets  because  this  material  probably  has  formed  at  temperatures  below 
400  degrees  Kelvin  (23:437).  There  is  some  evidence  that  some  asteroids 
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did  not  form  in  the  asteroid  belt,  but  at  much  greater  distances  from  the 
sun  as  comets,  [f  cometary  ice  is  trapped  in  the  interiors,  there  may  be 
an  outgassing  of  water  vapor  or  other  gasses  (23:437). 

The  only  data  that  is  now  available  on  the  chemical  composition  of 
asteroids  is  from  ground-based  observations.  By  assuming  that  the  sur¬ 
faces  of  the  asteroids  are  composed  of  common  cosmic  materials,  in  approx¬ 
imately  the  same  proportions,  spectral  studies  can  predict  the  likely 
composition  of  individual  asteroids  (16:25-40).  However,  dark  and/or 
opaque  material  remains  hidden  from  spectroscopic  study.  The  only  fea¬ 
sible  method  of  validating  these  remote  observations  is  to  compare  the 
results  to  the  composition  of  fragments  of  meteorites  found  on  the  Earth. 
Occasionally,  through  analyses  of  impact  craters  and  photographic  tracks 
of  meteorites,  it  is  possible  to  associate  them  with  an  asteroid  or 
comet.  Then  a  spectroscopic  comparison  can  be  made.  Just  as  likely,  a 
spectroscopic  comparison  can  infer  the  origin  of  a  meteorite  without  any 
o^her  supporting  evidence.  Unlike  the  Moon,  where  samples  returned  from 
the  surface  provide  positive  spectroscopic  proofing  of  remote  observa¬ 
tions,  asteroids  lack  any  confirmed  comparisons.  Therefore,  analysis  of 
surface  materials  remains  an  uncertain  proposition. 

The  exact  origin  of  meteorites  that  may  have  come  from  asteroids  is 
of  interest  for  other  reasons.  In  1890,  the  Farmington  meteroite  fell  in 
Kansas.  Although  no  photographic  evidence  was  available,  numerous  news¬ 
paper  and  eyewitness  reports  enabled  its  orbit  to  be  reconstructed. 
Analysis  revealed  that  the  meteorite  was  in  a  small  Earth  orbit  prior  to 
impact.  Further  back-plotting  of  the  orbit  associated  the  meteorite  with 
either  the  asteroid  Apollo,  Hermes  or  Cerebus  (15:234).  Direct  compari¬ 
son  of  the  L-chondrite  material  of  the  Farmington  meteorite  with  material 
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from  these  or  similar  asteroids  would  be  a  break-through.  Besides  the 
obvious  chemical-related  answers,  invaluable  information  could  be  obtained 
on  the  accuracy  and  reliability  of  some  of  the  orbital  paths,  and  the 
methods  for  establishing  the  origins  of  meteorites.  Some  researchers 
(22:54)  believe  that  most  of  the  meteorites  falling  on  the  Earth  are 
fragments  of  high-eccentricity  Apollo  asteroids  that  have  collided  with 
ether  objects  as  they  pass  through  the  heavily-populated  asteroid  belt. 
However,  no  direct  evidence  is  available  to  confirm  this  theory. 

The  most  likely  asteroids  to  have  originated  as  comets  are  those 
whose  orbits  have  a  high  eccentricity  and/or  a  high  inclination.  Hahn 
and  Rickman  (9:420)  tabulate  14  asteroids  (ranging  from  1036  Ganymed  to 
recently-discovered  1984  BC)  which  have  eccentricities  greater  than  0.47 
and  inclinations  up  to  34.57  degrees.  Cockrane  and  Barker  (4:296)  con¬ 
tend  that  most  of  the  numbered  asteroids  have  eccentricities  lower  than 
typical  comets,  and  have  diameters  that  are  too  large.  However,  the 
Aten,  Apollo  and  Amor  class  asteroids  are  in  orbits  much  more  closely 
matching  comets.  Many  of  these  asteroids  have  diameters  of  five  kilo¬ 
meters  or  less  in  diameter,  which  is  typical  of  comets.  Cockrane  and 
Barker  point  out  that  the  recently-discovered  asteroid  1983  TB  is  in  a 
"cometary"  orbit  that  is  almost  identical  to  the  orbit  of  19  Geminid 
meteors.  Spectral  studies  of  1983  TB  could  not  detect  any  gas  emission 
but  could  not  conclusively  deny  cometary  origin  (4:297).  Even  if  some 
asteroids  were  found  to  have  a  comet-like  composition,  there  remains  the 
orbital  dynamics  question  of  how  sufficient  numbers  of  comets  were  per¬ 
turbed  into  asteroid  orbits.  Apollo-Amor  asteroids  typically  have 
aphelion  distances  well  inside  the  orbit  of  Jupiter.  Periodic  comets 
typically  have  aphelion  distances  well  beyond  the  orbit  of  Jupiter. 


Yamada  (25:459)  hypothesizes  that  the  orbits  of  Apollo-Amor  asteroids 
could  be  perturbed  by  encounters  with  the  Earth  and  Venus  to  reduce  the 
large  aphelion  distances.  A  positive  link  between  comets  and  asteroids 
(such  as  a  material  sample)  would  promote  a  re-think  of  the  possible 
mechanisms  of  these  orbital  perturbations.  Farquhar  and  Dunham  (6:1) 
describe  realistic  trajectories  to  return  samples  from  Halley's  comet. 
While  it  is  too  late  for  these  missions,  they  do  portray  a  scenario  that 
could  be  used  with  other  comets.  The  lowest  cost  missions  (due  to  lowest 
energy  requirements)  are  fly-by  as  opposed  to  rendezvous.  A  proposed 
mission  would  collect  cometary  dust  during  the  high  speed  fly-by  and 
return  to  a  Shuttle-accessible  Earth  orbit  using  atmospheric  braking. 
While  this  technique  could  not  return  samples  from  an  asteroid,  it  does 
provide  for  part  of  the  puzzle  with  material  from  a  comet. 

Besides  the  scientific  interest  in  the  chemical  composition  of 
asteroids,  there  are  economic  interests.  If  just  one  asteroid  could  be 
explored,  the  baseline  data  could  be  used  to  determine  the  expected 
mineral  content  of  similar  asteroids.  Analysis  of  light  from  asteroids 
has  already  led  scientists  to  speculate  on  mineral  content  (23:434). 
Asteroids  with  similar  reflectivity  and  spectral  light  could  contain 
similar  minerals.  If  valuable  minerals  were  found,  it  might  be  more 
economical  (due  to  the  low  surface  gravity)  to  mine  an  asteroid  than  to 
mine  the  Moon.  If  strategical ly-important  minerals,  such  as  titanium  or 
uranium  were  found,  asteroid  exploration  could  have  significant  military 
impl ications . 

It  is  possible  that  the  orbit  of  a  small  asteroid  could  be  delib¬ 
erately  altered.  Since  many  of  the  recently-discovered  asteroids  are 
very  small,  often  less  than  500  meters  in  diameter  (11:44),  it  is 


conceivable  that  an  asteroid  could  be  manipulated  as  a  weapon.  It  may 
even  be  necessary  to  defend  against  a  random  asteroid  impact  (22:54). 
Wetherill  calculated  that  about  once  every  100  years,  an  asteroid  about 
one  kilometer  in  diameter  can  be  expected  to  pass  nearer  to  the  Earth 
than  the  distance  of  the  Moon.  Such  a  body  collides  with  the  Earth  about 
every  250,000  years.  The  energy  released  by  such  an  impact  would  equal 
the  yield  of  10,000  10-megaton  hydrogen  bombs.  However,  the  probability 
of  an  asteroid  actually  hitting  the  Earth  is  very  small  (8:171). 

Despite  the  remote  possibility  of  a  collision,  the  potential  consequences 
may  warrant  at  least  a  cursory  look  at  the  options.  French  and  Hulkower 
(8:167)  refer  to  three  methods  of  avoiding  a  collision.  The  orbit  of 
the  asteroid  could  be  altered  by  explosive  devices,  or  by  mass  drivers 
accelerating  asteroidal  material  to  produce  reaction  thrust.  For  a 
small  asteroid,  total  destruction  might  be  possible.  In  the  unlikely 
event  that  an  asteroid  were  detected  on  a  collision  course  with  Earth, 
current  knowledge  about  asteroids  would  make  it  very  difficult  to  mount 
one  of  these  defenses.  The  smaller  Earth-approaching  asteroids  with 
diameters  in  the  order  of  one  kilometer  are  the  least  understood,  but  of 
course  are  the  best  candidates  for  collisions.  No  accurate  size  or  mass 
data  is  available  so  the  magnitude  of  the  task  to  deflect  such  an  object 
is  based  on  inaccurate  information.  Composition  of  the  asteroid  would  be 
critical  if  the  mass  driver  or  total  destruction  options  were  contemplated. 
For  Earth-approaching  asteroids,  observational  data  to  determine  composi¬ 
tion  is  scant  or  nonexistent  (8:167).  Unfortunately,  if  information  and 
data  needed  to  deflect  an  asteroid  is  ever  required,  sufficient  time  to 
obtain  this  information  probably  will  not  be  available. 


Further  knowledge  of  the  size  and  mass  of  asteroids  is  being  gained 
by  observing  their  interaction  with  large  planets.  Williams  (24:1)  was 
able  to  achieve  significant  mass  determinations  of  the  three  largest 
asteroids  by  observing  the  perturbations  they  produced  in  the  orbit  of 
Mars.  About  three  dozen  other  asteroids  perturb  Mars  by  more  than  five 
meters  but  ranging  data  is  not  sufficiently  precise  to  produce  accurate 
mass  determinations.  Ranging  data  available  from  Viking  spacecraft  can 
be  used  to  determine  the  motion  of  the  inner  planets  against  an  inertial 
frame  of  reference.  A  major  objective  of  this  ranging  is  to  check  on  the 
constancy  of  the  gravitational  constant.  However,  for  further  Mars 
missions  to  contribute  to  this  check,  more  accurate  information  on  the 
sizes  and  densities  of  asteroids  will  be  necessary  (24:7).  The  accuracies 
required  can  only  be  determined  accurately  by  direct  measurement  during 
fly-by  or,  preferably,  rendezvous  missions. 

With  the  renewed  interest  in  the  Earth-approaching  asteroids,  there 
has  been  a  corresponding  increased  interest  in  the  orbital  mechanics 
required  to  plan  and  support  any  possible  rendezvous  mission.  Most  of 
the  capability  of  conducting  meaningful  analyses  of  intercept  trajectories 
and  long  term  projections  has  resided  on  large  mainframe  computers.  With 
the  advent  of  powerful  microcomputers,  and  the  development  of  new  algo¬ 
rithms  to  analyze  trajectories,  accurate  and  timely  analyses  are  now 
possible  using  modest  computing  resources. 


III.  Objectives 


The  overall  objective  of  this  paper  is  to  demonstrate  a  computer¬ 
ized  analysis  of  the  accessibility  of  Earth-approaching  asteroids  with 
particular  emphasis  on  the  three  most  accessible  asteroids,  and  those 
that  have  been  discovered  from  January  to  November  1985. 

To  accomplish  the  overall  objective,  the  sub-objectives  are: 

1.  To  determine  which  asteroids  have  been  discovered  during  the 
period  of  interest.  (Earth-approaching  asteroids  have 
recently  been  discovered  at  the  rate  of  between  three  and 
eight  per  year.  As  well,  previous  discoveries  are  occasionally 
added  or  deleted  from  the  list  of  these  asteroids  as  the 
orbital  parameters  are  refined  by  subsequent  observations.) 

2.  To  accurately  calculate  the  orbits  of  the  Earth  and  the 
selected  asteroids,  and  the  intercept  trajectories. 

3.  To  determine  the  accessibility  of  these  asteroids. 

4.  To  compare  the  accessibility  of  these  asteroids  with  previous 
discoveries. 

5.  To  compare  the  capability  of  the  computer  programs  with  results 
generated  by  main  frame  computer  facilities. 

6.  To  extend  the  current  knowledge  of  the  accessibility  of  Earth- 
approaching  asteroids. 


Ciassi fication  of  Asteroids 

To  be  designated  as  Earth-approaching,  asteroids  obviously  must 
approach  the  Earth.  The  Apollo  class  asteroids  are  in  orbits  very 
similar  to  the  orbit  of  the  Earth,  and  occasionally,  the  orbital  paths 
come  close  to  intersecting.  A  new  and  significant  asteroid  discovery  was 
made  on  a  photograph  taken  on  the  night  of  27/28  February  1982  at  the 
Palomar  Mountain  Observatory.  The  asteroid,  designated  as  1982  DB, 
passed  within  4.08  million  kilometers  of  the  Earth  about  one  month 
before  its  discovery.  An  analysis  of  the  orbit  of  1982  DB  revealed  that 
it  was  an  Apollo  asteroid  (an  asteroid  whose  orbit  crosses  the  Earth's 
orbit).  This  discovery  was  particularly  significant  because  the  asteroid 
had  passed  nearer  the  Earth  than  any  other  asteroid  since  1976,  and  its 
orbit  made  it  the  easiest  of  all  known  asteroids  to  reach  with  a  probe. 

A  detailed  analysis  of  intercept  opportunities  from  1982  to  2002 
revealed  that  not  only  is  1982  DB  easy  to  reach,  but  that  two-way  mission 
opportunities  occur  almost  yearly  (11:42). 

Besides  the  Apollo  asteroids  which  cross  the  Earth's  orbit,  there 
are  two  other  classes  of  asteroids  that  occasionally  approach  the  Earth. 
Aten  asteroids  orbit  the  sun  just  inside  the  Earth's  orbit,  and  Amor 
asteroids  orbit  just  outside.  All  three  classes  could  include  asteroids 
that  would  provide  favorable  rendezvous  opportunities.  Collectively, 
these  asteroids  are  known  as  AAA  ( Apol lo-Aten-Amor)  or  Earth-approaching 
asteroids. 

“An  asteroid  is  said  to  be  Earth-approaching  if  its  distance 

at  perihelion,  q,  is  within  1.3  AU  of  the  sun,  thereby  admitting 


the  possibility  of  passing  within  a  few  tenths  of  an  All  of  Earth. 
These  bodies  are  divided  into  three  families  according  to  q  and 
semimajor  axis,  a: 


Family 

Orbital  Characteristics 

Aten 

q  <=  1.02, 

a  <=  1 

Apollo 

q  <=  1.02, 

a  >  1 

Amor 

1.02  <=  q 

<=  1.3  " 

The  "family"  names  are  derived  from  a  representative  asteroid  in 
each  family  whose  orbit  satisfies  these  conditions.  Nominal  orbital 
parameters  of  the  AAA  namesakes  for  epoch  12h  00  Dec  1985  f 14:3)  are: 


Number  Name 

q 

a 

e 

i 

AN 

AP 

M 

2026  ATEN 

.79 

.97 

.18 

19 

108 

148 

12 

1862  APOLLO 

.65 

1.47 

.56 

6 

35 

285 

326 

1221  AMOR 

1.09 

1.92 

.43 

12 

171 

26 

50 

The  Palomar  Planet-Crossing  Asteroid  Survey  has  discovered  several 
of  these  Earth-approaching  asteroids  since  1982.  However,  many  of  the 
orbital  parameters  are  still  only  approximations,  and  only  a  few  of  these 
new  asteroids  have  yet  been  given  an  official  catalog  number.  Until  they 
have  been  observed  on  several  different  orbits,  and  their  orbital  para¬ 
meters  accurately  determined,  they  will  be  known  only  by  the  year  of 
discovery  followed  oy  a  two-letter  serial  code.  Despite  this  uncer¬ 
tainty,  the  initial  orbital  parameters  are  generally  sufficiently  accu¬ 
rate  to  calculate  how  easily  and  how  often  a  rendezvous  could  be  made. 

Optimal  Trajectories 

Mathematical  techniques  have  been  devised  to  calculate  the  energy 
required  to  rendezvous  with  a  given  asteroid  at  a  given  time.  However, 
at  a  different  time,  the  intercept  geometry  can  be  vastly  aifferent.  It 
may  not  be  at  all  apparent  whether  or  not  this  new  geometry  represents  a 
better  or  worse  intercept  opportunity.  Each  instant  in  time  can  be 


analyzed  directly  and  accurately,  but  no  direct  method  has  yet  been 
devised  to  find  the  optimum  intercept  opportunities.  It  is  currently 
believed  that  no  set  of  equations  can  be  found  that  will  directly  solve 
this  proDlem,  and  no  attempt  will  be  made  in  this  paper  to  do  so. 

Techniques  have  recently  been  developed  to  indirectly  solve  the 
minimum  energy  intercept  problem  (17:1-5;  12:458-461).  These  techniques, 
practical  only  when  implemented  on  a  computer,  require  that  intercept 
energy  be  calculated  for  every  conceivable  position  of  an  asteroid  in 
its  orbit  for  every  conceivable  position  of  the  Earth  in  its  orbit.  The 
total  number  of  calculations  required  depends  on  the  desired  resolution 
in  the  answer.  If  the  positions  of  each  body  were  taken  each  degree 
around  its  360  degree  orbit,  a  total  of  129,600  (360  x  360)  calculations 
would  be  required  for  each  asteroid.  Each  of  these  calculations  is  far 
from  trivial.  For  every  intercept  opportunity,  the  orbital  position  of 
the  Earth  and  the  asteroid  must  be  calculated,  as  well  as  the  satellite 
trajectory  near  the  Earth,  and  the  interplanetary  trajectory  under  the 
influence  of  the  sun.  Parts  of  these  calculations  cannot  be  solved  in 
closed  form,  out  require  computer  iterative  techniques  to  get  an  approxi¬ 
mate  answer.  Because  of  the  large  number  of  calculations  involved,  even 
a  computer  can  spend  long  hours  trying  to  find  the  answers. 

Fortunately,  for  the  purpose  of  the  analyses  conducted  in  this  paper, 
only  a  few  of  the  3500  known  asteroids  are  classified  as  Earth-approaching 
In  1983,  McFadden  listed  36  Amors,  34  Apollos  and  3  Atens  (16:25).  Many 
of  these  73  asteroids  are  obviously  less  accessible  than  two  that  have 
received  recent  attention.  Detailed  calculations  have  been  done  for 
1943  Anteros,  and  1982  DB.  The  latter  was  found  to  be  the  most  accessible 
Since  McFadden's  list  in  1983,  there  have  been  nine  new  Earth-approaching 


asteroids  discovered,  and  two  or  three  more  are  being  found  each  year. 


Accessibi 1 ity 

The  accessibility  of  an  asteroid  is  the  quantitative  assessment  of 
the  ease  by  which  an  asteroid  may  be  approached  or  explored.  Ease  is 
quantified  by  a  combination  of  distance  from  the  Earth,  orbital  orien¬ 
tation,  orbital  position  and  velocity,  and  gravitational  force.  The 
further  an  asteroid  is  from  the  Earth,  the  greater  the  energy  and  time 
required  to  reach  it.  An  orbit  that  is  at  a  highly-inclined  angle  to  the 
Earth's  orbit  may  be  several  times  more  difficult  to  reach  than  an  orbit 
that  is  in  the  same  plane.  The  relative  positions  and  velocity  differ¬ 
ences  between  the  two  orbits  determine  when,  and  how  often,  travel  is 
practical  from  a  body  in  one  to  a  body  in  the  other.  The  gravitational 
force  of  an  asteroid  determines  the  amount  of  energy  that  must  be  expended 
to  arrive  and  depart  from  the  surface  of  the  asteroid.  In  summary,  many 
factors  must  be  considered  to  determine  the  accessibility  of  an  asteroid. 

Because  of  the  complexity  of  the  accessibility  calculations,  only 
the  most  promising  of  intercept  opportunities  have  been  analyzed  in 
detail,  and  then  only  for  a  period  of  about  20  years.  Detailed  analyses 
has  been  done  for  some  of  the  recently-discovered  Earth-approaching 
asteroids.  These  analyses  need  to  be  extended  far  enough  into  the  future 
to  detect  any  unusually  favorable  opportunities.  This  time  period  should 
include  any  opportunity,  or  situation  (such  as  possible  defensive  action 
or  an  ambitious  rendezvous  mission)  for  which  planning  would  have  to 
be  started  in  the  next  several  years.  Typically,  three  (occasionally  as 
many  as  eight)  new  Earth-approaching  asteroids  are  discovered  each  year. 
These  discoveries  regularly  provide  new  candidates  for  these  analyses. 


Methodology 

The  accessibility  of  selected  asteroids  was  determined  by  using  an 
assimilation  of  classical  orbital  mechanics  and  recently-developed  mathe¬ 
matical  techniques.  Because  of  the  complexity  of  the  calculations,  com¬ 
puter  programming  was  used  extensively  throughout  this  project.  Calcu¬ 
lations  and  computer  programs  were  validated  by  comparing  the  results 
to  known  results.  No  simulated  or  arbitrary  data  was  used.  All  input 
data  was  real  and  the  most  recent  that  is  available.  All  output  was  cal¬ 
culated  as  accurately  as  practical  for  the  techniques  used.  Details  of 
computer  accuracy  are  listed  in  the  section  on  computer  facilities. 

There  are  numerous  techniques  for  the  calculation  of  rendezvous 
trajectories  and  the  accessibility  of  asteroids.  Both  closed-form 
equations  and  iterative  methods  are  available.  Most  of  these  approaches 
provide  the  accuracy  but  not  the  speed  to  handle  the  large  number  of 
calculations  required  to  determine  accessibility.  Main  belt  asteroids 
(all  those  in  similar  near-circular  orbits  between  Mars  and  Jupiter)  can 
be  assessed  with  a  few  typical  calculations  that  will  not  vary  signifi¬ 
cantly  from  one  to  tne  other.  However,  Earth-approaching  asteroids  each 
present  unique  opportunities  with  highly-varying  accessibilities.  More 
efficient  calculations  are  needed. 

Two  new  techniques  were  announced  by  Ross  and  Hulkower  in  1981,  and 
were  demonstrated  in  the  calculation  of  the  accessibility  of  asteroid 
Anteros  (17:1-5).  The  first  method  determined  the  optimal  rendezvous 
trajectories  and  optimal  time  of  flight  for  each  relative  position  of 
the  Earth  and  the  asteroid.  The  second  technique  used  a  plot  of  these 
optimal  trajectories  for  all  combinations  of  positions.  These  plots  were 


called  Prime  Rib  curves  because  of  the  contour  patterns  produced.  A 
further  routine  was  developed  to  increase  the  accuracy  and  efficiency  of 
finding  the  best  of  all  of  these  trajectories.  In  1983,  Ross  and 
Hulkower  made  a  more  detailed  analysis  of  missions  (including  return 
missions)  to  Anteros.  Later  in  1983,  an  asteroid  which  had  been  dis¬ 
covered  in  February  1982  (named  1982  DB) ,  was  reported  to  have  replaced 
Anteros  as  the  most  accessible  known  asteroid  (11:42-46).  For  Anteros 
and  1982  DB,  the  published  analyses  list  a  total  of  36  rendezvous  tra¬ 
jectories,  and  show  representative  Prime  Rib  curves  against  which  some 
of  the  calculations  in  this  paper  are  compared. 

Initially,  the  Ross  and  Hulkower  technique  appeared  to  be  the  most 
promising.  However,  extensive  analysis  of  the  fundamental  equations  and 
hundreds  of  runs  of  computer  programs  revealed  a  number  of  problems  with 
this  approach.  The  primary  problem  was  a  consistent  reluctance  of  the 
iterative  equations  to  converge.  Techniques  to  enable  convergence  for 
one  set  of  data  proved  to  be  unstable  for  new  data.  When  the  equations 
converged,  reasonable  results  were  obtained,  but  when  the  equations  did 
not  converge,  no  results  were  obtained.  Because  the  overall  analysis 
depended  on  1296  independent  runs  of  the  algorithms  to  produce  one  plot 
of  accessibility  for  one  asteroid,  it  was  vital  that  the  computer  be 
able  to  run  automatically  with  no  concern  about  equations  converging.  It 
wa s  then  learned  from  Hulkower  that  he  too  had  trouble  with  these  equa¬ 
tions.  In  a  recent  paper  by  Hulkower,  Lau  and  Bender  (12:458-461),  he 
abandoned  the  original  approach,  and  offered  a  completely  new  algorithm. 
The  Prime  Rib  curves  were  retained  as  the  best  presentation  of  the 
optimal  trajectories. 
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The  calculations  of  the  accessibility  of  asteroids  are  based  on  the 


Hulkower,  Lau  and  Bender  techniques.  Orbits  of  the  Earth  and  selected 
asteroids  were  calculated  by  computer  for  the  period  5  January  1985  to 
5  January  2020.  These  orbits  provide  the  framework  in  which  to  apply 
the  analytical  techniques.  Primary  attention  was  to  be  given  to  rendez¬ 
vous  missions,  but  two-way  and  fly-by  missions  were  also  addressed. 
Particularly  interesting  trends  beyond  1  January  2010  were  sought.  The 
output  of  optimal  trajectories  and  launch  opportunities  was  compared 
with  results  calculated  by  Ross  and  Hulkower,  Hulkower,  Lau  and  Bender, 
and  with  other  sources  (18:1-16).  Calculations  for  asteroids,  scenarios 
and  times  not  available  from  other  sources  are  presented  in  the  appendices 
of  this  paper. 

The  following  steps  summarize  the  development  of  the  algorithms: 

1.  A  PASCAL  computer  program  was  developed  to  calculate  the  orbits 
of  the  Earth  and  the  asteroids  from  05  January  1985  to 

02  January  2020. 

2.  The  orbital  program  was  validated  against  ephemeris  data  from 
Astronomical  Almanac  for  the  Year  1985  and  Almanac  for 
Computers  1985. 

3.  The  orbital  program  was  extended  to  produce  positions  and 
velocities  of  the  Earth  and  the  asteroids  in  a  form  usable  by 
the  subsequent  optimal  trajectory  algorithms. 

4.  Initially,  the  analytical  techniques  of  Hulkower  and  Ross 
were  converted  to  PASCAL  subroutines. 

5.  Orbital  parameters  of  the  asteroid  Anteros  were  used  to 
exercise  the  subroutines. 

6. '  Prime  Rib  curves  for  Anteros  were  compared  to  results  pub¬ 

lished  by  Hulkower  and  Ross. 

7.  The  Hulkower  and  Ross  algorithms  were  aDandoned  in  favor  of 
the  new  ones  developed  by  Hulkower,  Lau  and  Bender. 

8.  New  PASCAL  computer  subroutines  were  written.  The  original 
orbital  programs  and  the  Prime  Rib  plotting  routines  were 
reusable. 


9.  The  new  algorithms  were  shown  to  be  stable,  accurate  and 
efficient. 

10.  A  PASCAL  program  was  developed  to  associate  Julian  Dates  and 

true  anomalies  of  the  Earth  and  the  asteroids  to  2  January  2020. 
This  facilitated  finding  optimal  launch  opportunities  from 
Prime  Rib  plots. 

The  Ross  and  Hul kower  Technique 

As  an  extension  to  Battin's  work  on  optimal  one  impulse  transfers 
between  circular  orbits,  Ross  and  Hulkower  developed  a  technique  to  solve 
one  and  two  impulse  transfers  between  orbits  of  arbitrary  eccentricity 
(17:1-5).  For  ease  of  reference,  this  technique  will  be  referred  to  as 
the  Ross  technique.  Ross  used  the  same  geometric  representation  as 
Sattin  and  others  as  depicted  in  Figure  1. 


Figure  1.  Geometry  of  the  Rendezvous  Problem 


Position  vectors  r^  and  locate  the  departure  planet  and  arrival 

planet  in  relation  to  the  sun.  These  vectors  are  joined  by  a  vector  c 

representing  the  vector  subtraction  of  r^  from  r^.  The  angle  theta  is 

measured  from  r^  to  r^,  positive  in  the  direction  of  rotation  of  r^ ,  and 

ranging  from  0  to  360  degrees.  Ross  establishes  three  unit  vectors  for 

r^,  r^  and  c,  and  then  handles  all  his  equations  in  a  new  reference  frame 

defined  by  sums  and  differences  of  these  three  unit  vectors.  Departing 

slightly  from  his  notation,  the  unit  vectors  are  labelled  U  ,  Ur  and  U 

rl  '  2  C 

Positions  and  velocities  are  then  calculated  as  multiples  of  the  vectors 
(U  +  Ur  ),  (U  -  U  )  or  (U  -  Up  ),  (U  +  U_  ).  Although  this  trans- 

L  I  ^  U  l  2  *2 

formation  simplifies  the  equations  used  to  solve  the  optimal  impulse 
problem,  it  considerably  increases  the  manipulation  required  to  handle 
these  equations  on  a  computer.  The  following  equations  necessary  to 
implement  the  computer  solution  are  extracted  from  the  much  more  detailed 


development  by  Ross  (17:2-3),  and  Battin  (2:93-112): 


Vj  =  K[q  sece(gc  +  U^)  +  SG  tan  e  (Uc  -  U  )] 


V*  =  a*(Uc  +  U  )  +  B*(UC  -  U  )  +  Y*  h 


V,  =  -K[q  sec  e(U_  -  U  )  +  SG  tan  e(Ur  +  U_  )] 

L  2 


\  *  "p<uc  -  Ur  >  *  s.<uc  *  ur  >  *  vp  ** 


k  =  /uc/[2s(s-c)] 


(1) 


(2) 


(3) 


(4) 


(5) 


c 


Pv  " 


SG  =  +1  for  transfer  orbits  not  passing  through  apoapsis 
between  7^  and  TV, 

SG  =  -1  for  transfer  orbits  passing  through  apoapsis  between 
7^  and  Tv, 

e  =  free  variable  in  iterations 
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The  final  equation  (11)  was  derived  from  previous  Ross  equations  but 
differs  significantly  from  his  version  of  the  same  result.  It  appears 
that  there  is  at  least  a  typographical  error,  and  probably  a  development 
error  in  the  final  Ross  equation.  The  differences  include  the  quantity 
2k  and  several  plus  and  minus  signs.  His  equation  does  not  collapse 
back  to  his  previous  result  for  the  one-impulse  case  as  it  should  by 
equating  k2  to  zero.  As  well,  his  equation  would  not  converge  to  produce 
an  optimal  impulse  answer.  When  it  was  learned  that  a  number  of  problems 
had  plagued  this  algorithm,  and  that  a  new  one  had  been  developed  by  one 
of  the  original  authors,  this  approach  was  abandoned. 


The  Hulkower,  Lau  and  Bender  Technique 


Hulkower,  Lau  and  Bender  developed  a  patched  conic  iterative  method 
of  determining  the  optimum  two-impulse  transfers  between  arbitrary 
elliptical  orbits  (12:458-461).  For  ease  of  reference,  this  technique 
will  be  referred  to  as  the  Hulkower  technique.  Although  this  technique 
is  not  as  elegant  as  the  Ross  technique,  it  is  somewhat  simpler  to  imple¬ 
ment  on  a  computer.  All  positions  and  velocities  are  calculated  in 
heliocentric-ecliptic  coordinates  thus  requiring  no  transformations  into 
unusual  frames  of  reference.  As  with  the  Ross  method,  the  transfer  begins 
in  a  circular  low  Earth  orbit  and  ends  in  an  elliptical  or  circular  orbit 
around  the  asteroid.  Variations  of  this  basic  transfer  are  allowed  for 
fly-by  and  for  return  missions. 

Hulkower  begins  by  defining  the  changes  of  velocity  at  the  departure 
and  arrival  points  in  terms  of  the  hyperbolic  excess  velocity,  the  gravi¬ 
tational  constants  of  the  departure  and  arrival  planets,  and  the  periapse 
and  apoapse  d: stances  of  the  departure  and  arrival  parking  orbits. 


=  [T \  +  (2u1/q1)]i  -  {2u1Q1/Cq1(Q1  +  q1)]}i  (12) 

dV2  =  CT|  +  (2u2/q2)]i  -  {2p2Q2/[q2(Q2  +■  q2)]}J  (13) 


whe-e 

A V  =  velocity  change 
T  =  hyperbolic  excess  velocity 
u  =  gravitational  constants  of  planets 
q  =  perigee  distance 
Q  =  apogee  distance 


Hulkower  adds  the  two  changes  in  velocity  to  produce  one  expression 
for  the  total  change  of  velocity.  He  then  takes  the  partial  derivative 
of  this  expression  with  respect  to  the  parameter  p  of  the  transfer  orbit, 
producing  an  expression  that  equals  zero  at  the  minimum  total  change  of 
velocity.  Although  this  expression  is  not  difficult  to  handle  in  the 
computer  program,  it  was  not  used  in  favor  of  an  iterative  form  of  the 
derivative.  However,  it  can  be  used  to  confirm  that  the  optimal  trajec¬ 
tory  has  been  found.  Ross  does  not  detail  his  method  of  finding  the 
optimal  trajectories  but  it  probably  used  this  expression. 

Following  previous  work  by  McCue  and  Bender  (12:459),  Hulkower 
writes  an  expression  relating  velocities  on  the  transfer  orbit  to  the 
hyperbolic  excess  velocity,  the  velocity  of  the  departure  or  arrival 
planet,  unit  position  vectors,  and  a  defined  velocity  v.  Separate  equa¬ 
tions  are  necessary,  depending  on  the  size  of  the  angle  between  the  posi¬ 
tion  vectors  of  the  departure  and  arrival  bodies.  Following  Hulkower, 
the  upper  sign  refers  to  angles  less  than  180  degrees,  and  the  lower 
sign  refers  to  angles  between  180  and  360  degrees. 


±  ( v  +  z  Uj ) 

(14) 

i  (V  +  z  U2) 

(15) 

This  defined  velocity  is  expressed  in  terms  of  the  gravitational 
constant  of  the  sun,  the  parameter  of  the  transfer  orbit,  and  the  posi¬ 
tion  vectors  of  the  departure  and  arrival  planets  in  the  heliocentric- 
ecliptic  reference  frame. 


V  =  [(up)*  (K  -  ?,)/  j 7.  x  "r  | 


(16) 


The  angle  between  the  position  vectors  is  used  to  calculate  a  quan¬ 


tity  z. 

z  =  (u/p)  *  tan($/2)  (17 

where 

0  =  angle  between  "r,  and 

The  *ree  variable  that  is  used  to  find  the  optimal  transfer  trajec¬ 
tory  is  the  parameter  of  the  transfer  orbit.  Hulkower  gives  expressions 
for  the  upper  and  lower  bounds  of  the  parameter  to  simplify  the  search 
process.  These  maximum  and  minimum  values  of  the  parameter  are  given  in 
terms  of  the  position  vectors  of  the  departure  and  arrival  points. 


rlr2  -  rl  '  r2 


min 


rl  *  r2  *  C2(rlr2  *  rl  '  r2>] 
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rlr2  rl  '  r2 


"max 


rrv  [2(rir2  *  ri  •  r2)] 
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Hulkower  plots  the  optimal  changes  in  velocity  using  the  same  Prime 
Rib  curves  that  were  used  for  the  Ross  method.  Each  point  shows  the 
total  change  of  velocity  for  a  given  combination  of  true  anomaly  of  the 
departure  planet  at  launch  and  the  true  anomaly  of  the  arrival  planet  at 
arrival.  Figure  2  shows  the  essential  elements  of  a  Prime  Rib  curve. 
Figure  3  shows  an  example  of  the  Prime  Rib  plot  generated  by  the  program 
in  this  paper.  This  particular  plot  is  for  the  asteroid  1982  DB,  the 
most  accessible  asteroid  yet  discovered. 


ASTEROID  1982DB 


V .  Computer  Program  Devel opment 
Common  PASCAL  Procedures 

GetOate.  This  procedure  is  used  to  obtain  both  the  start  and  the 
end  date,  and  the  time  of  an  ephemeris  period.  The  input  must  be  properly 
entered  or  the  procedure  will  reject  it  and  ask  again.  The  month  must 
contain  three  letters,  with  the  first  letter  upper  case  and  the  remain¬ 
ing  two  letters  lower  case.  The  day  of  the  month  is  error  checked  and 
must  be  an  integer  between  one  and  the  maximum  days  for  that  month.  The 
year  must  be  an  integer  between  1900  and  2099,  the  limits  of  the  Julian 
date  algorithm  used  in  the  program  (20:B2).  The  hour  of  the  day  is  the 
Coordinated  Universal  Time  (UTC)  input  as  a  real  number  between  0.0  and 
24.0.  The  procedure  returns  the  year,  month,  day  and  hour  as  real 
numbers . 

CalendarOate.  This  procedure  accepts  a  Julian  date  and  returns  a 
real  number  representation  of  the  day,  month  and  year.  Leap  years  are 
handled  correctly. 

Jul lanDate.  This  procedure  is  the  complement  to  Cal endarDate.  It 
accepts  a  year,  month,  day  and  decimal  hour  (UTC)  and  returns  a  real 
number  Julian  Date.  The  full  seven  digits  plus  the  decimal  part  of  the 
Julian  date  are  returned.  The  algorithm  is  from  the  Almanac  for  Computers 
1985  (20:B2) . 

Sol veKeplersEqn .  This  procedure  accepts  a  mean  anomaly  in  degrees 
and  returns  the  eccentric  anomaly  in  degrees.  It  is  implemented  using  a 
traditional  iterative  approach,  with  a  selectable  accuracy. 


SunXYZCoord .  This  procedure  accepts  the  Julian  date  and  uses  the 
global  variables  of  the  Earth's  orbital  elements  in  the  heliocentric- 
ecliptic  coordinate  system.  Outputs  are  the  X,Y,Z  coordinates  of  the  sun 
in  the  geocentric-equatorial  coordinate  system.  The  obliquity  of  the 
ecliptic  is  a  variable  adjusted  from  a  reference  date.  This  algorithm 
was  adapted  from  the  low-precision  formulae  for  the  Sun's  coordinates 
from  the  Astronomical  Almanac  1985  ( 20 : C24 ) .  The  stated  precision  is 
0.01  degree  between  the  years  1950  and  2050.  This  formula  sets  the 
limits  over  wnich  this  program  can  be  used.  Greater  accuracy  would 
require  either  vastly  more  complex  formulae,  or  a  series  expansion  tech¬ 
nique  such  as  that  used  to  generate  ephemerides  in  the  Astronomical 
Almanac.  The  latter  is  very  accurate  but  can  only  be  used  over  the 
period  for  which  the  coefficients  for  the  expansion  are  tabulated  in  the 
almanac  (usually  one  year). 

GetDataFi le.  This  procedure  handles  the  data  file  manipulation  for 
most  of  the  programs  used  in  this  paper.  The  data  files  are  stored  as 
text  files  so  they  can  easily  be  viewed,  corrected  or  updated.  All  the 
data  files  are  formatted  in  the  same  way.  They  contain  the  name  and/or 
number  of  the  solar  system  object  plus  the  classical  orbital  element  set 
referenced  to  the  heliocentric-ecliptic  coordinate  system.  Although  the 
programs  which  use  this  procedure  will  normally  be  concerned  with  the 
Earth  and  Earth-approaching  asteroids,  this  procedure  was  written  to 
accept  a  variety  of  data  files.  All  nine  planets  are  included  in  the 
file  PLANET. ELE,  which  contains  the  1985  orbital  elements.  Rendezvous 
missions  could  be  calculated  for  Mars  and  Venus  to  compare  the  trajec¬ 
tories  and  the  energy  requirements  with  those  for  missions  to  the 
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asteroids.  Several  well-known  comets  are  included  in  a  file  called 
COMET. ELE.  The  main  purpose  of  this  file  was  to  verify  the  ephemesis 
program  by  comparing  the  output  to  published  ephemerides  of  Halley's 
comet.  The  similarity  of  some  cometary  orbits  to  those  of  asteroids 
permit  comparisons  of  rendezvous  trajectories.  Several  files  are  avail¬ 
able  for  the  low-numbered  asteroids,  again  because  these  objects  are 
well  observed  and  ephemerides  are  available.  Specific  comparisons  were 
made  of  published  and  generated  ephemerides  for  asteroids  1  Ceres, 

2  Pallas,  3  Juno  and  4  Vesta. 

The  input  menu  for  this  procedure  suggests  using  one  of  the  standard 
files  for  planets,  comets  or  asteroids,  and  allows  an  optional  user 
defined  file  to  be  used.  For  some  programs  used  in  this  paper,  the  user 
defined  option  was  replaced  by  a  specific  file  of  AAA  asteroids. 

InitialDisplay.  This  procedure  generates  the  initial  display  for 
the  main  program.  It  introduces  the  program  and  lists  which  data  files 
are  available  on  the  data  disk.  It  also  controls  the  input  of  the  name 
of  the  optional  user-supplied  data  file  if  that  option  was  selected.  No 
error  checking  is  done  to  see  if  a  user-supplied  data  file  exists  before 
access  is  attempted,  so  the  user  must  ensure  that  the  proper  data  file 
is  on  the  default  disk  drive. 


PASCAL  Program  Hel i oCentri cQrbi t 

The  program  Hel ioCentricOrbit  is  used  to  generate  ephemerides  for 


solar  system  objects  (1:51-83;  7:90-98).  The  program  computes  right 


and  assemble  the  data  files.  Any  or  the  data  files  can  be  used,  permit¬ 
ting  the  ephemerides  of  planets,  comets  and  asteroids  to  be  generated. 

The  period  of  validity  of  the  program  is  1950  to  2050.  The  accuracy 
depends  on  the  algorithm  for  the  position  of  the  sun  discussed  in  the 
procedure  SunXYZCoord,  and  the  effects  of  the  perturbations. 

CalculateObjPosition.  This  procedure  is  the  heart  of  the  program. 

It  first  determines  the  position  of  the  sun  in  the  geocentric-equatorial 
reference  frame  by  calling  the  procedure  SunXYZCoord.  The  classical 
element  set  is  then  used  to  calculate  the  position  of  the  solar  system 
object  in  the  heliocentric-ecliptic  reference  frame.  Traditional  methods 
are  used  in  these  calculations,  such  as  the  use  of  auxiliary  quantities 
and  Gaussian  constants.  The  procedure  to  solve  Kepler's  equation  is 
called  to  calculate  the  eccentric  anomaly.  All  positions  are  transformed 
into  the  geocentric-equatorial  frame  and  the  right  ascension  and  decli¬ 
nation  are  generated  from  tngometric  calculations. 

PeriodEphemeris.  This  procedure  permits  the  user  to  establish  the 
period  for  which  the  ephemeris  is  to  be  calculated,  or  optionally,  a 
single  date  and  time.  If  a  period  of  time  is  selected,  the  user  is  asked 
for  the  start  and  end  dates  and  times,  and  the  interval  of  time  between 
each  calculation.  All  the  dates  are  converted  to  seven  digit  Julian 
dates  for  ease  of  calculations.  Output  may  be  directed  to  either  the 
screen  or  printer.  The  same  format  is  output  to  whichever  display  is 
selected.  The  output  consists  of  the  object  name,  the  day,  month  and 
year,  and  the  right  ascension  and  declination  in  hours  and  degrees.  This 
procedure  uses  two  internal  procedures  called  CalcEphemeris  and  Print  to 
control  the  calculations  and  to  generate  the  output. 


PASCAL  Program  RossHul kowerDel taVee 

This  program  uses  the  Ross  and  Huikower  algorithms  to  generate  a 
plot  of  optimal  trajectory  change  in  velocity  or  delta  vee.  For  sim¬ 
plicity,  this  technique  will  be  referred  to  as  the  Ross  technique.  Parts 
of  the  program  are  derived  from  the  program  Hel ioCentricOrbit.  This  pro¬ 
gram  is  self  contained,  except  for  the  data  files.  It  generates  posi¬ 
tions  and  velocities  for  both  the  departure  and  arrival  planets  and  runs 
the  optimal  trajectory  algorithms.  Output  can  be  either  a  Prime  Rib  plot 
of  delta  vee  for  each  combination  of  true  anomalies  of  the  departure  and 
arrival  planets,  or  a  tabular  listing  of  the  same  information. 

Several  functions  such  as  "tan"  or  tangent  are  programmed  in  this 
program  because  they  are  not  part  of  the  fundamental  PASCAL  compiler. 
These  functions  are  derived  from  the  sine,  cosine  and  arctangent  func¬ 
tions  using  standard  trigometric  relations.  No  further  explanations  will 
be  given  for  these  and  similar  functions  as  their  implementation  is  self- 
evident  in  the  program  listing. 

Because  the  procedures  and  algorithms  in  this  program  must  be 
repeated  1296  times  to  generate  one  Prime  Rib  plot,  efficiency  is  an 
important  consideration.  Several  of  the  procedures  from  the  program 
Hel ioCentricOrbit  have  been  modified  to  increase  their  efficiency. 
Quantities  that  have  to  be  computed  only  once  and  remain  constant  for  a 
number  of  iterations  are  often  broken  out  into  separate  procedures. 
Modifications  of  the  previous  procedures  to  accomplish  this  efficiency 
are  obvious  in  the  program  listings  and  will  not  be  discussed  further. 

CaiculateObj Posit ion.  This  procedure  is  identical  to  the  procedure 
of  the  same  name  in  the  program  Hel ioCentricOrbit. 


Cal culateObj Velocity.  This  procedure  calculates  the  vector  veloci¬ 
ties  of  the  departure  and  arrival  planets.  These  vector  quantities  are 
computed  from  the  positions,  orbital  elements  and  heliocentric  gravita¬ 
tional  constant  using  standard  formulae  (13:300). 

Departure360.  This  procedure  uses  the  previous  two  procedures  to 
calculate  the  position  and  velocity  vectors  for  the  departure  planet 
every  ten  degrees  of  true  anomaly  from  0  to  360  degrees.  The  ten  degree 
increment  is  arbitrary  and  is  selected  only  because  it  is  convenient. 

This  selection  provides  a  compromise  between  speed  of  execution  and 
resolution  of  output.  Any  increment  could  be  substituted  without  affect¬ 
ing  the  proper  functioning  of  the  program,  but  a  smaller  quantity  (such 
as  one  degree)  could  greatly  increase  the  run  time  of  the  program.  If 
the  output  is  to  the  Prime  Rib  plot,  provision  must  be  made  for  the 
different  resolution.  A  practical  compromise  would  be  to  run  the  pro¬ 
gram  first  with  ten  degree  resolution  to  determine  areas  of  interest,  and 
rerun  the  program  in  higher  resolution  over  a  limited  combination  of  true 
anomalies.  The  output  of  this  procedure  is  stored  in  two  2-dimensional 
arrays  of  position  and  velocity  for  the  departure  planet. 

Arrival360.  This  procedure  functions  exactly  like  the  previous  one. 
The  reason  for  having  a  separate  procedure  for  the  departure  and  arrival 
planets  is  to  facilitate  operating  the  program  over  a  limited  range  of 
true  anomalies  that  probably  will  be  different  for  each  planet. 

GetAl phaBetaGamma .  This  procedure  converts  the  position  and  velocity 
data  from  the  heliocentric-ecliptic  reference  frame  into  a  skewed  frame 
used  by  the  Ross  algorithms.  The  Alpha,  Beta  and  Gamma  are  used  to  be 
consistent  with  the  Ross  notation,  and  refer  to  the  three  components  of 


velocity.  In  the  arrays,  [1],  [2],  and  [3]  correspond  to  Alpha,  Beta 
and  Gamma. 

The  transformations  are  computed  by  the  dot  product  of  the  velocity 
vectors  and  a  unit  vector  in  the  new  frame  as  established  by  Ross.  Once 
the  transformation  has  been  completed,  the  calculations  are  predomi¬ 
nantly  2-dimensional  in  the  plane  defined  by  the  transfer  orbit.  Three- 
dimensions  are  used  when  the  results  are  changed  back  to  the 
heliocentric-ecliptic  frame.  The  angular  momentum  of  the  departure 
planet  and  the  cross  product  of  the  two  position  vectors  are  used  to 
establish  the  positive  orientation  for  the  orbit  direction. 

UnitVectorcrlr2.  This  procedure  is  key  to  the  transformation 
required  in  procedure  GetAl phaBetaGamma .  It  provides  the  unit  vectors 
that  define  the  new  reference  frame  established  by  Ross  (17:2).  This  new 
basis  is  orthogonal,  but  the  values  of  Alpha,  Beta  and  Gamma  are  not 
direct  transormations  to  these  unit  vectors.  Rather,  they  are  referred 
to  the  sum  and  difference  of  pairs  of  these  unit  vectors.  This  causes 
the  quantities  Alpha,  Beta  and  Gamma  to  be  distorted  versions  of  the 
true  velocity  vectors.  In  other  words,  the  magnitude  of  the  velocity 
after  the  transformation  is  not  the  same  as  the  magnitude  of  the 
velocity  before  the  transformation.  This  must  be  taken  into  considera¬ 
tion  when  these  quantities  are  converted  back  to  the  original  reference 
frame. 

GetNextSetTrueAnomal ies.  This  procedure  accesses  the  arrays  con¬ 
taining  the  positions  and  velocities  for  the  departure  and  arrival 
planets  in  preparation  for  calculating  the  Alpha,  Beta  and  Gamma  pro¬ 
gram,  and  accessed  as  needed  by  the  delta  vee  part  of  the  program.  This 


permits  more  flexibility  in  limiting  the  combinations  of  true  anomalies 
over  which  the  delta  vee  algorithms  will  run.  As  well,  it  permits  the 
program  to  be  more  modular  in  design  to  facilitate  modifications. 

FindDeltaV.  This  subroutine  is  the  heart  of  the  Ross  algorithm  to 
find  the  optimal  trajectories.  All  the  previous  procedures  were  merely 
setting  the  stage  and  assembling  the  data  for  this  procedure.  All 
variables  required  by  this  procedure  are  passed  to  it  by  the  calling 
statement.  No  global  variables  are  accessed  from  within  the  procedure. 
This  permits  FindDeltaV  to  operate  as  a  unit  independent  from  the  rest 
of  the  program.  As  such,  it  could  be  easily  transported  to  other  pro¬ 
grams.  The  main  part  of  this  procedure  begins  by  defining  several  con¬ 
stants  according  to  Ross.  The  entire  procedure  is  a  double  iteration. 

The  solution  ultimately  requires  that  the  quantities  DeltaVl  and  DeltaV2 
be  found.  These  quantities  are  used  to  calculate  k,  and  which  are  in 
turn  used  to  calculate  epsilon.  Epsilon  is  then  used  to  calculate 
DeltaVl  and  DeltaV2,  and  the  cycle  repeats  until  the  desired  accuracy  is 
reached.  The  equations  for  k^  and  k^  incorporate  the  size  of  the  parking 
orbit  and  the  gravitational  constant  for  both  the  departure  and  arrival 
planets.  This  permits  these  orbits  to  be  arbitrarily  chosen.  As  well, 
by  equating  k^  to  zero,  no  orbit  is  established  at  the  arrival  planet, 
and  the  result  is  a  fly-by  mission. 

FindEpsi Ion .  This  procedure  is  used  to  calculate  epsilon  for  a 
given  ^  and  k,,.  Because  epsilon  cannot  be  isolated  on  one  side  of  the 
equation,  an  iterative  solution  is  required.  The  iteration  is  very 
sensitive  to  initial  conditions,  and  can  diverge  and  produce  no  solution. 
When  it  does  converge,  it  does  so  very  rapidly.  Some  of  the  problem 


with  the  converging  can  probably  be  attributed  to  the  differences  in  the 
Ross  formula  for  epsilon  and  the  one  derived  in  this  paper.  The  Ross 
formula  appears  to  have  at  least  typographical  errors.  The  quantity 
(2k)  is  missing  from  the  denominators,  and  several  of  the  plus  and  minus 
signs  are  switched.  Deriving  this  equation  from  Ross'  previous  equations 
produced  the  equation  used  in  this  procedure.  Ross'  equation  for  Epsilon 
would  not  function  at  all. 

Delta360.  This  procedure  controls  the  procedures  that  arrange  the 
data  and  calculate  the  delta  vee.  It  effectively  separates  the  first 
part  of  the  program  which  calculates  the  positions  and  velocities  from 
the  second  part  which  implements  Ross'  algorithms. 

PASCAL  Program  Hul kowerLauBenderDel taVee 


This  program  was  developed  to  replace  the  RossHul kowerDel taVee,  and 
was  used  to  generate  all  the  optimal  trajectory  Prime  Rib  plots,  tabular 
form  of  the  Prime  Rib  plots,  and  the  times-of-fl ight  associated  with  each 
trajectory.  The  program  uses  the  Hulkower,  Lau  and  Bender  algorithms 
(12:458-461).  For  simplicity,  Hulkower  will  be  used  to  mean  Hulkower, 

Lau  and  Bender.  These  algorithms  were  found  to  be  highly-stable  and 
accurate  for  all  combinations  of  true  anomalies,  and  for  all  the  combi¬ 
nations  of  departure  and  arrival  planets  that  were  tried.  The  algorithms 
presented  by  Hulkower  were  changed  and  simplified  for  the  computer  program 
with  no  apparent  loss  of  stability  or  accuracy.  These  changes  are  docu¬ 
mented  in  the  explanations  of  the  individual  procedures.  As  with  the 
previous  programs,  use  was  made  of  the  common  procedures  for  file  input 
and  initial  calculations  to  set  up  the  data  for  the  Hulkower  algorithms. 
Some  changes  were  made  to  these  procedures  to  account  for  slightly 


different  data  requirements,  but  these  changes  are  minor,  and  are  readily 
apparent  in  the  PASCAL  code.  No  additional  explanation  will  be  given 
here. 

Pi ff I .  This  procedure  is  the  heart  of  the  implementation  of 
Hulkower' s  algorithms.  It  uses  three  internal  procedures  called  FindDiff, 
Newton  and  GaussGodalTimeOfFl ight  which  are  unique  to  this  program.  The 
procedure  accesses  the  global  variables  for  departure  and  arrival  planet 
position  and  velocity  vectors,  and  transfers  the  data  to  local  variables 
using  the  Hulkower  notation.  Several  quantities  that  remain  constant  are 
precomputed  prior  to  calling  the  other  procedures.  The  parameter  of  the 
transfer  orbit  is  used  as  the  free  variable  that  is  iterated  to  find  the 
optimal  trajectory.  To  reduce  the  number  of  iterations,  and  to  provide 
stability  to  the  search  routines,  upper  and  lower  bounds  are  computed  for 
the  parameter.  The  starting  point  for  the  search  routines  then  becomes 
the  value  of  the  parameter  halfway  between  its  computed  maximum  and  mini¬ 
mum  value.  The  Hulkower  algorithms  identify  two  cases  of  the  transfer 
orbit  geometry  that  must  be  handled  differently.  The  cases  are  deter¬ 
mined  by  the  size  of  the  angle  between  the  position  vectors  of  the  depar¬ 
ture  and  arrival  planets.  If  the  angle  is  less  than  or  equal  to  180 
degrees,  then  one  set  of  equations  apply.  If  the  angle  is  greater  than 
180  degrees,  then  a  second  set  of  equations  apply.  Fortunately,  the  only 
difference  between  the  two  sets  of  equations  is  the  plus  or  minus  signs 
for  a  number  of  terms  in  the  equations.  For  the  computer  implementation, 
a  quantity  call  SL  is  set  to  either  +1.0  or  -1.0  and  is  inserted  through¬ 
out  the  equation.  The  optimal  trajectory  is  computed  for  the  appropriate 
case.  This  feature  of  the  algorithm  provides  a  check  on  its  proper 
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functioning.  If  both  cases  are  run  through  the  optimal  trajectory  pro¬ 
cedure,  the  case  which  gives  the  minimum  delta  vee  should  be  the  one 
that  is  appropriate  to  the  size  of  the  angle  between  the  position  vectors 
of  the  departure  and  arrival  planets.  This  check  for  the  proper  func¬ 
tioning  of  the  algorithms  was  retained  in  the  program  that  is  listed  in 
Appendix  A.  However,  it  is  not  essential  to  the  operation  of  the  pro¬ 
gram  and  could  be  removed  to  provide  an  increase  in  speed  of  the  program. 

Newton .  This  procedure  is  called  when  the  preparatory  calculations 
described  in  the  previous  procedure  have  been  completed.  The  purpose  of 
Newton  is  to  iterate  the  parameter  of  the  transfer  orbit,  beginning  at  a 
point  half  way  between  the  calculated  extremes  and  remaining  at  all  times 
between  these  extremes.  For  a  given  value  of  the  parameter,  the  procedure 
FindDiff  is  called  to  find  the  difference  between  the  delta  vee  for  that 
value  of  the  parameter  and  the  delta  vee  for  the  Drevious  iteration  of 
the  parameter.  The  value  of  the  parameter  for  the  next  iteration  is 
determined  by  halving  the  previous  iteration  step  and  moving  either  plus 
or  minus  depending  on  the  result  of  the  previous  two  iterations.  When 
the  difference  in  the  total  delta  vee  converges  to  a  defined  accuracy, 
the  iteration  terminates  and  the  delta  vee  for  that  value  of  the  parameter 
is  noted.  To  accommodate  the  situation  where  the  iteration  oscillates 
between  two  very  small  values,  a  limit  of  50  iterations  is  maintained. 
Convergence  normally  occurs  after  about  30  iterations,  and  is  not  very 
sensitive  in  accuracy.  Near  the  optimal  point  of  the  curve  of  this 
function,  the  slope  is  very  flat  and  gradual,  and  high  accuracies  are 
easy  to  attain. 


GaussGodalTimeOfFl ight.  This  procedure  is  designed  to  take  the  few 
bits  of  information  about  the  transfer  orbit  that  have  so  far  been  found 
and  calculate  the  time  of  flight  of  the  transfer  trajectory.  Gauss  and 
Godal  each  developed  formulae  to  solve  the  time  of  flight  problem  (2:82- 
84).  However,  each  solution  depended  on  knowing  information  about  the 
orbit  that  is  not  available  from  the  Hulkower  algorithms.  By  combining 
parts  of  the  Gauss  and  Godal  solutions,  an  equation  was  found  that  included 
only  known  terms.  Although  this  equation  is  a  large  and  awkward  one,  it 
is  effective  and  stable  for  all  orbits  investigated  in  this  paper.  Care 
must  be  taken  with  the  computer  implementation  of  this  equation  to  prevent 
overflew  and  excessive  round-off  errors  because  it  contains  several  cubes 
and  fourth-roots  of  transcendental  expressions.  Inputs  to  the  procedure 
are  four  scaler  quantities  (the  magnitude  of  the  two  position  vectors 
of  the  departure  and  arrival  planets,  the  angle  between  these  vectors, 
and  the  parameter  of  the  transfer  orbit).  The  output  is  the  time-of- 
flight.  With  the  time-of-fl ight  known,  all  the  other  elements  of  the 
transfer  orbit  can  be  calculated  if  required.  Therefore,  this  procedure 
provides  the  important  link  between  delta  vee  plots  of  the  Hulkower 
algorithms  and  complete  knowledge  of  the  transfer  orbit. 

PASCAL  Program  CalcAnomal ies 

This  program  is  designed  to  calculate  the  true  anomalies  of  both 
the  departure  and  arrival  planets  every  ten  days  from  05  January  1985  to 
02  January  2020.  It  generates  a  data  file  called  AN0M-1.DTA  that  can  be 
accessed  by  the  program  ListAnomal ies.  This  program  does  not  incorporate 
a  procedure  to  access  the  data  files  containing  the  orbital  elements  of 
all  the  selected  solar  system  objects.  The  element  set  for  the  Earth  is 


included  but  the  element  set  for  the  target  must  be  provided.  The  pro¬ 


gram  was  designed  like  this  so  it  could  readily  be  used  as  a  procedure 
within  an  all-encompassing  program.  Besides  generating  the  data  file, 
this  program  simultaneously  lists  the  true  anomalies  in  the  same  format 
as  the  program  List Anomalies.  The  reason  that  two  separate  programs  were 
written  is  that  the  data  files  can  be  accessed  much  faster  than  the  data 
can  be  calculated.  It  is  therefore  prudent  to  calculate  this  data  once 
and  store  the  results.  This  is  particularly  important  when  the  data 
must  be  accessed  numerous  times  by  various  search  routines. 

KeplerEquation.  This  procedure  solves  the  classic  Keplerian  equa¬ 
tion  for  the  eccentric  anomaly  given  the  mean  anomaly.  The  algorithm 
was  adapted  from  one  presented  by  Danby  and  Burkardt  (5:95-107)  that 
features  quintic  convergence.  It  is  extremely  fast,  achieving  conver¬ 
gence  to  at  least  ten  decimal  places  in  three  iterations  or  less  in 
almost  all  cases.  It  continues  to  perform  at  this  level  even  at  eccen¬ 
tricities  greater  than  0.9999.  This  speed  ensures  that  this  procedure 
and  program  operate  efficiently  in  calculating  large  numbers  of  true 
anomal ies . 

Bacon .  This  procedure  calculates  the  true  anomalies  from  the 
element  set  using  the  procedure  KeplerEquation. 

SaveAnomalyFile.  This  procedure  generates  the  data  file  AN0M-1.DTA 
from  the  two  arrays  of  true  anomalies.  The  data  file  is  formatted  to  be 
compatible  with  the  program  Li stAnomal ies . 


PASCAL  Program  ListAnomal ies 

This  program  is  designed  to  read  a  data  file  AN0M-1.DTA  that  was 
generated  by  the  program  CalcAnomal ies.  It  will  list  true  anomalies  for 


both  the  departure  and  arrival  planet  every  ten  days  from  05  January 
1985  to  02  January  2020.  The  data  is  loaded  into  two  1-dimensional 
arrays  that  each  hold  1278  entries.  If  a  higher  resolution  of  days  was 
required  (such  as  every  two  days),  the  array  size  and  the  maximum  number 
of  anomalies  would  have  to  be  adjusted  to  match  those  established  by  the 
program  Cal cAnomal ies .  The  output  is  formatted  with  Julian  dates  in 
increments  of  100  days  in  the  left  column,  and  in  increments  of  10  days 
across  the  top. 

PASCAL  Program  FindLaunchDate 

This  program  finds  a  launch  date  (if  one  exists)  given  a  true  anomaly 
of  the  departure  planet  at  launch,  the  true  anomaly  of  the  arrival  planet 
at  arrival,  and  the  time-of-fl ight.  It  loads  a  data  file  created  by  the 
program  CalcAnomal ies  containing  the  true  anomalies  of  the  departure 
planet  and  the  arrival  planet  every  ten  days  from  January  1985  to 
January  2020.  The  search  routine  asks  for  the  two  true  anomalies  and 
the  time-of-fl ight.  It  then  expands  the  search  to  include  the  original 
true  anomalies  plus  and  minus  ten  degrees,  and  the  original  time-of- 
fl  ight  plus  and  minus  ten  days.  All  combinations  of  three  values  of 
three  quantities  are  included.  Thus,  for  each  set  of  inputs,  27 
searches  are  conducted  to  include  all  adjacent  values.  This  expansion 
is  done  to  account  for  the  resolution  of  ten  degrees  and  ten  days  with 
the  associated  round-off  errors.  It  also  permits  an  area  search  rather 
than  a  point  search.  The  validity  of  an  area  search  is  readily  seen 
from  the  Prime  Rib  plot  where  it  is  apparent  that  adjacent  values  of 
delta  vees  are  usually  quite  similar  for  values  of  true  anomalies  that 
vary  ten  degrees  or  more.  Output  from  this  program  converts  Julian 


dates  to  day,  month  and  year  (the  launch  date),  and  prints  the  associated 
true  anomalies  and  time-of-fl i ght  found  in  the  search  routine. 


LoadAnomalyFi le.  This  procedure  loads  the  1278  sets  of  true 
anomalies  for  ten  day  increments  from  January  1985  to  January  2020.  This 
data  is  found  in  the  data  file  created  by  the  program  CalcAnomal ies 
called  AN0M-1.DTA  or  in  a  user-supplied  file  of  the  same  format.  The 
true  anomalies  are  stored  in  the  file  as  real  numbers  with  decimals. 
LoadAnomalyFi le  converts  these  anomalies  to  integer  quantities  rounded 
to  the  nearest  ten  degrees.  The  integer  days  are  also  rounded  to  the 
nearest  ten  days. 

JtoD  and  MriteDate.  These  procedures  convert  the  Julian  date 
derived  from  the  array  position  and  a  constant  to  a  calendar  date.  The 
result  is  printed  when  needed  in  the  form  DD  Month  YYYY  (i.e.,  21  January 
1998) . 


InputTAandTQF.  This  procedure  asks  the  user  for  the  true  anomalies 
of  the  departure  and  arrival  planets  (in  degrees)  and  for  the  time-of- 

flight  in  days.  This  information  is  derived  from  the  Prime  Rib  and  the 

time-of-fl i ght  tables.  Typically,  the  user  would  select  a  combination 
of  true  anomalies  that  represented  a  minimum  value  of  delta  vee  on  the 

Prime  Rib  plot  and  the  program  would  respond  with  the  dates  if  any  when 

these  combinations  of  true  anomalies  and  time-of-fl ight  would  occur 
during  the  period  from  January  1985  to  January  2020. 

DoCombos.  This  procedure  and  its  embedded  procedures  expand  the 
search  area  to  include  three  values  of  each  of  the  three  inputs.  All  27 
combinations  are  then  checked  for  occurrence  during  the  time  period. 

All  matches  that  are  found  are  printed,  and  the  program  asks  if  further 
runs  are  required. 
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Computer  Resources 

These  programs  were  designed  to  operate  efficiently  on  a  modest 
microcomputer  system  using  a  widely-available  computer  language. 

The  programs  were  developed  on  an  IBM  Personal  computer  with  640 
kilobytes  of  internal  memory.  The  system  was  equipped  with  two  floppy 
disk  drives  each  with  a  capacity  of  360  kilobytes  per  disk.  An  8087  math 
co-processor  integrated  circuit  chip  was  installed  to  increase  the  speed 
of  mathematical  operations.  A  monochrome  monitor  displayed  25  lines  of 
30  characters.  The  monitor  was  also  used  as  a  monochrome  graphics 
terminal  with  a  resolution  of  640x200  pixels.  Printout  was  directed  to 
an  Epson  FX-80  dot  matrix  printer.  It  was  also  used  as  a  graphics 
terminal  to  reproduce  the  640x200  plots  displayed  on  the  monochrome 
monitor. 

All  programming  was  done  in  PASCAL  using  the  Borland  International 
Turbo  PASCAL  compiler/editor  system.  To  take  advantage  of  the  math  co¬ 
processor  chip.  Version  3.0  of  Turbo  PASCAL  with  8087  support  was  used 
for  all  programs. 

A  minimum  system  to  support  these  programs  would  require  256  kilo¬ 


bytes  of  internal  memory,  one  disk  drive  and  a  monochrome  display.  The 
programs  are  usable  (with  some  inconvenience)  without  the  graphics  out¬ 
put,  without  the  speed  of  the  math  co-processor  chip,  and  without  the 


VI.  Results  and  Discussion 


The  results  of  this  analysis  of  the  accessibility  of  Earth- 
approaching  asteroids  consist  of  computer  programs  to  implement  the  algo¬ 
rithms,  Prime  Rib  plots  and  tables  to  show  the  global  optimum  Delta  Vee 
information  for  selected  asteroids,  and  tables  of  launch  opportunities 
for  selected  asteroids  from  the  year  1985  to  2020. 

Listings  of  the  PASCAL  code  for  the  computer  programs  are  given  in 
Appendices  A,  B,  C,  D,  E  and  F.  These  listings  are  important  because 
they  represent  a  very  practical  tool  for  anyone  with  an  interest  in  the 
accessibility  of  not  only  the  asteroids  analyzed  in  this  paper,  but  for 
the  analysis  of  any  asteroid,  comet  or  planet  in  the  solar  system.  The 
orbital  elements  of  selected  solar  system  objects  are  given  in  Appendix 
M.  The  volume  of  output  chat  would  be  generated  to  analyze  just  a  few 
asteroids  is  very  large.  Only  a  selected  sample  of  the  more  interesting 
and  recently-discovered  asteroids  can  be  presented  here.  Therefore,  the 
computer  programs  are  a  must  to  do  a  detailed  analysis  of  a  particular 
target  of  interest.  These  programs  also  provide  a  tool  to  quickly 
analyze  future  discoveries.  Since  these  discoveries  have  been  occurring 
at  the  rate  of  between  three  and  eight  per  year  in  recent  years,  and  are 
likely  to  continue  with  the  renewed  interest  in  this  area,  these  computer 
programs  will  provide  a  readi ly-accessible  source  of  detailed  study  as 
soon  as  the  initial  orbital  elements  of  the  discoveries  are  known. 

The  Prime  Rib  plots  and  tables  provide  a  detailed  look  at  the  global 
picture  of  the  accessibility  of  selected  asteroids  in  a  form  that  is 
easy  to  comprehend.  They  show  not  only  the  specific  information  for  one 


asteroid  but  can  be  viewed  together  to  quickly  appreciate  the  relative 
accessibility  of  a  number  of  asteroids.  These  Prime  Rib  plots  are  some¬ 
what  different  from  the  Prime  Rib  curves  presented  by  Ross  and  Hulkower 
in  that  they  employ  shading  rather  than  contour  lines  to  convey  the 
information.  The  shading  has  the  advantage  that  it  is  readily  apparent 
which  plots  or  segments  of  plots  portray  the  best  potential  opportunities 
Very  simply,  the  darker  the  plot,  the  lower  the  required  delta  vee.  The 
plots  with  the  largest  areas  of  dark  shading  further  indicate  the  aster¬ 
oids  for  which  the  most  launch  opportunities  are  likely  to  exist.  The 
plots  are  much  easier  to  generate  than  the  curves  and  can  readily  be 
produced  on  a  dot-matrix  printer  that  is  a  part  of  most  modest  computer 
systems . 

The  tables  of  launch  opportunities  cover  the  period  from  1985  to 
2020  for  the  three  most  accessible  asteroids.  These  tables  retain  the 
true  anomaly  increment  at  10  degrees  while  Lau  and  Hulkower  (14:1-27) 
refined  the  increment  considerably,  producing  additional  launch  oppor¬ 
tunities.  This  paper  introduces  a  detailed  analysis  of  three  new  Earth- 
approaching  asteroids  that  have  been  discovered  in  1985.  Launch  oppor¬ 
tunities  for  these  new  asteroids  are  also  tabulated  from  1985  to  2020. 


Computer  Programs 

This  section  is  a  summary  of  the  function  of  computer  programs  that 
have  been  developed  to  support  the  analyses  presented  in  this  paper. 
Further  details  will  be  apparent  in  the  easily-read  PASCAL  code  in 
Appendices  A  through  F.  Each  program  is  designed  to  operate  interac¬ 
tively  with  the  user  and  no  further  explanation  should  be  required  to 
use  each  effectively.  Suggested  extensions  to  the  programs  are  proposed 


In  the  following  chapter  of  this  paper.  Each  program  is  listed  by  the 
descriptive  name  that  is  embedded  in  the  PASCAL  code. 

Program  Function  and  Brief  Description 

1.  Hel ioCentricOrbit  This  is  a  validation  program  that 

demonstrates  the  accuracy  of  the  algo¬ 
rithms  that  calculate  the  orbital 
positions  of  solar  system  bodies. 

This  program  is  embedded  in  both  the 
Delta  Vee  programs  to  provide  the 
fundamental  position  and  velocity 
data  required  by  the  Delta  Vee  algo¬ 
rithms.  A  sample  of  the  accuracy  of 
this  program  compared  to  published 
Astronomical  Almanac  ephemerides  is 
tabulated  in  Appendix  l. 

2.  RossHul kowerDel taVee  This  program  implements  the  algo¬ 

rithms  introduced  by  Ross  and 
Hulkower  (17:1-5).  As  explained 
previously,  these  algorithms  were 
found  to  be  unstable  and  difficult  to 
handle.  This  program  is  presented 
because  of  the  suggestions  offered 
in  the  next  chapter  of  this  paper. 


3.  Hul kowerLauBenderDel taVee  This  program  is  by  far  the  most  impor¬ 

tant  to  the  analyses  discussed  in 


this  paper.  It  produces  the  Prime 
Rib  plots  and  tables  that  depict  the 
global  accessibility  of  asteroids. 

It  also  produces  a  corresponding 
time-of-fl ight  for  each  trajectory 
shown  in  the  Prime  Rib  plots.  The 
output  from  this  program  is  used  by 
the  following  programs  to  generate 
launch  opportunities. 

The  strong  points  of  this  pro¬ 
gram  are  its  stability,  its  accuracy, 
its  speed  and  its  flexibility. 
Throughout  its  many  runs  on  a  variety 
of  asteroids  and  other  bodies,  it  was 
100%  stable  and  showed  none  of  the 
problems  of  iterations  failing  to 
converge  that  were  found  in  the 
RossHulkowerDeltaVee  program.  Com¬ 
parisons  with  Hulkower,  Lau  and 
Bender  data  showed  an  accuracy  of 
about  0.01  kilometers  per  second  for 


delta  vees  of  between  5.0  and  60.0 


kilometers  per  second,  and  accuracies 
of  between  1  and  20  days  for  times- 


of-f light  up  to  700  days.  Suggested 
improvements  to  these  accuracies  are 
discussed  in  the  following  chapter. 


4.  CalcAnomal ies 


5.  Li stAnomal ies 


This  program  operates  very  quickly 
considering  the  massive  amounts  of 
iterations  required  by  the  algorithms. 
One  Prime  Rib  plot,  a  table  of  delta 
vees  and  table  of  times-of-fl ight 
showing  1296  individual  rendezvous 
trajectories  are  produced  in  a  total 
time  of  15.5  minutes.  This  time  is 
constant  for  any  asteroid.  The  flexi¬ 
bility  of  this  program  permits  it  to 
use  any  planet,  comet  or  asteroid  in 
orbit  around  the  sun  as  either  depar¬ 
ture  or  arrival  body.  Rendezvous, 
fly-by,  and  return  missions  are  all 
supported. 


This  program  calculates  true  anomalies 
for  a  pair  of  bodies  (usually  the 
Earth  and  an  asteroid)  every  ten  days 
from  5  January  1985  to  2  January 
2020.  It  operates  in  less  than  a 
minute  producing  a  file  that  is  used 
by  the  following  programs. 


This  program  produces  a  listing  of 
the  data  generated  by  CalcAnomal ies. 
It  operates  at  the  speed  of  the 
printer  used  for  output. 
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6.  FindLaunchOate  This  program  finds  any  launch  oppor¬ 

tunities  that  exist  for  a  given  com¬ 
bination  of  true  anomalies  for  the 
period  5  January  1985  to  2  January 
2020.  The  accuracy  within  this  pro¬ 
gram  is  plus  or  minus  five  days  based 
on  one  half  the  increment  of  the 
search.  Despite  having  to  investi¬ 
gate  27  combinations  of  true  anomalies 
and  times-of-fl ight  for  a  period  of 
35  years,  it  operates  with  no  percep- 
table  delay. 

Prime  Rib  Plots  and  Tables 

Prime  Rib  plots  are  displayed  for  asteroids  1982DB,  1982XB,  Anteros, 
1985JA,  1985PA  and  1985TB  in  Appendix  G.  Tabulated  values  of  the  plots 
and  times-of-fl ight  are  very  lengthy.  Therefore,  only  the  values  for 
one  selected  asteroid  are  shown.  The  asteroid  selected  is  1985TB  because 
it  was  the  most  recent  and  the  most  accessible  of  the  three  Earth- 
approaching  asteroids  discovered  in  1985.  The  values  for  1985TB  are 
detailed  in  Appendices  H  and  I.  The  following  table  summarizes,  in 
order  of  accessibility,  the  global  optimum  trajectory  information  con¬ 
tained  i n  the  detailed  tables  for  all  six  asteroids.  The  table  also 
shows  the  relative  accessibility  of  these  asteroids  compared  to  the  76 
that  have  been  ranked. 


1982XB 

Anteros 


1985JA 

1985PA 


TABLE  I 

Accessibility  of  Selected  Asteroids 


Global 

Minimum 

True  Anomaly 

True  Anomaly 

Total 

of  Earth 

of  Asteroid 

Time- 

Delta  V 

at  Launch 

at  Arrival 

of-Fl ight 

Asteroid  (km/s) 

(deg) 

(deg) 

(days)  Rank 

20 

150 

229 

1 

120 

185 

2 

140 

240 

435 

3 

320 

110 

218 

61 

120 

260 

407 

66 

60 

220 

302 

73 

One  of  the  common  characteristics  of  all  of  the  Prime  Rib  plots  are 
the  vivid  lines  of  demarcation  separating  the  plots  into  three  areas. 
Since  the  plots  wrap  around  both  horizontally  and  vertically,  the  lines 
actually  separate  the  plots  into  two  areas  separated  by  180  degrees. 
These  areas  show  how  the  accessibility  changes  throughout  the  360  degree 
range  of  true  anomalies.  If  only  the  true  anomalies  of  the  departure 
planet  are  considered,  one  would  expect  that  there  would  be  points  on 
this  range  where  the  accessibility  is  best  and  points  where  it  is  worst. 
The  lines  of  demarcation  are  the  worst  areas  and  between  these  lines 
there  are  a  total  of  two  best  points,  one  usually  being  significantly 
better  than  the  other.  Like  the  worst  areas,  the  best  areas  are  sepa¬ 
rated  by  about  180  degrees.  As  the  phase  of  the  target  true  anomaly 
changes  with  respect  to  the  departure  true  anomaly,  the  lines  of  demarca 
tion  propagate  across  the  plot  at  an  angle  that  averages  45  degrees. 

For  complex  geometries,  such  as  one  orbit  being  highly  inclined  to  the 


other,  or  one  orbit  being  much  larger  than  the  other,  the  plots  become 
more  complex  with  more  local  minimum  and  maximum  areas. 

Prime  Rib  Plot  for  1982DB.  This  plot  shows  the  global  optimal  point 
at  about  20  degrees  of  true  anomaly  of  the  Earth  and  150  degrees  for  the 
asteroid.  A  second  local  minimum  occurs  at  20  degrees  and  220  degrees. 
Although  both  areas  have  about  the  same  value  of  delta  vee,  the  first 
area  is  considerably  larger  than  the  second.  Size  of  the  area  has  a 
direct  impact  on  the  size  of  the  launch  window.  The  larger  the  area  for 
a  given  delta  vee,  the  larger  will  be  the  launch  window  to  take  advantage 
of  that  particular  delta  vee.  Just  as  important,  the  larger  the  area, 
the  higher  the  probability  that  the  Earth  and  the  asteroid  will  actually 
find  themselves  at  that  combination  of  true  anomalies.  If  the  area  is 
very  small,  such  as  a  degree  or  so  across  (and  beyond  the  resolution  of 
this  particular  plot),  it  will  probably  be  tens  or  hundreds  of  years 
before  that  particular  combination  of  true  anomalies  will  occur.  If 
the  asteroid  were  in  an  orbit  with  a  period  that  was  an  exact  multiple 
of  the  period  of  the  Earth,  then  only  a  few  of  the  combinations  of  true 
anomalies  would  ever  occur.  This  would  mean  that  parts  of  the  Prime  Rib 
plots,  although  possibly  interesting,  would  be  of  no  value  for  mission 
planning.  If  the  global  optimal  point  occurred  in  this  region,  it  would 
never  be  achievable  because  the  asteroid  and  the  Earth  simply  would 
never  get  into  the  right  positions.  For  1982DB,  this  problem  does  not 
exist.  The  global  optimal  condition,  or  a  point  very  close  to  it  occurs 
frequently.  Almost  the  entire  plot  for  1982DB  shows  delta  vee  values  of 
less  than  20  kilometers  per  second,  and  large  areas  less  than  7  kilometers 
per  second.  These  large  areas  of  small  delta  vees  would  indicate  that 


1982DB  is  a  very  accessible  asteroid.  In  fact,  combined  with  reasonable 
titnes-of-fl ight  and  f  equent  launch  windows,  these  values  of  delta  vees 
make  19820B  the  most  accessible  asteroid  yet  discovered. 

Prime  Rib  Plot  for  1982XB.  This  plot  is  remarkably  similar  to  the 
one  for  1982DB,  except  that  it  is  displaced  slightly  along  the  axis  of 
the  true  anomaly  of  the  Earth.  The  global  optimal  delta  vee  of  about 
5.2  kilometers  per  second  makes  1982XB  about  0.7  kilometers  per  second 
less  accessible  than  1982DB,  but  still  the  second  most  accessible 
asteroid . 

Prime  Rib  Plot  for  Anteros.  The  plot  for  Anteros  shows  global 
optimal  values  just  slightly  worse  than  1982XB.  Like  the  other  two, 
Anteros  displays  large  areas  of  relatively  low  delta  vee  enabling  the 
combinations  of  true  anomalies  which  represent  good  launch  opportunities 
to  occur  frequently.  Anteros  is  ranked  third  in  accessibility. 

Prime  Rib  Plot  for  1985JA.  This  plot  immediately  shows  that  1985JA 
is  much  less  accessible  than  the  three  most  accessible  asteroids.  Only 
the  lightest  shade  of  grey  was  plotted,  indicating  no  values  of  delta 
vee  below  15  kilometers  per  second.  Most  of  the  plot  is  white  meaning 
that  these  values  of  delta  vee  are  too  high  to  be  of  any  real  interest. 

If  these  higher  values  were  of  interest,  the  pi f  .  could  be  replotted 
with  the  lowest  delta  vee  normalized  to  black,  and  the  shades  of  grey 
extending  to  higher  values.  The  global  optimal  point  on  the  plot  is  15.7 
kilometers  per  second  giving  1985JA  a  ranking  of  66  for  accessibility  out 
of  76  that  are  ranked. 

Prime  Rib  Plot  for  1985PA.  No  values  of  delta  vee  plotted  for 
1985PA  with  the  default  setting  for  shades  of  grey.  The  global  optimal 


point  was  21.5  kilometers  per  second  ranking  the  asteroid  73rd  out  of  76, 
just  about  the  most  inaccessible  of  the  known  Earth-approaching  aster¬ 
oids.  This  large  value  is  due  mainly  to  the  nigh  inclination  of  55 
degrees  of  the  orbit.  Plane  changes  in  intercept  trajectories  are  very 
costly  in  delta  vee.  All  three  of  the  1985  asteroids  have  fairly  high 
inclinations  and  fairly  low  accessibility. 

Prime  Rib  Plot  and  Table  for  1985TB.  This  asteroid  is  the  most 
accessible  of  the  three  discovered  in  1985.  It  shows  a  plot  with  two 
shades  of  grey,  and  a  global  optimal  point  of  13.3  kilometers  per  second. 
That  still  places  it  only  61st  out  of  76  that  are  ranked.  A  fairly  small 
area  of  the  plot  shows  values  near  13.3,  so  launch  dates  near  that  value 
would  be  infrequent. 

Overall,  the  three  1985  asteroids  are  very  inaccessible  compared 
with  the  three  most  accessible.  They  require  about  three  to  five  times 
the  delta  vee  to  reach  than  the  most  accessible,  1982DB.  For  comparison, 
they  have  roughly  the  same  accessibility  as  the  main  belt  asteroid 
Pallas.  1985TB  is  slightly  more  accessible  than  Pallas,  1985JA  is  about 
the  same,  and  1985PA  is  slightly  worse.  That  comparison  demonstrates 
that  asteroids  that  approach  the  Earth  are  not  necessarily  more  accessible 
than  distant  asteroids.  Inclination  is  an  important  factor,  as  is  the 
shape  of  the  orbit  and  the  velocity  of  the  asteroid  when  it  is  near  the 
Earth. 

Launch  Opportunities  to  the  Three  Most  Accessible  Asteroids 


a  selection  of  the  optimum  and  near-optimum  trajectories  taken  from  the 
Prime  Rib  plots.  They  are,  of  course,  not  the  only  launch  opportunities 
because  given  enough  delta  vee,  one  can  launch  any  time.  Only  rendezvous 
missions  are  considered  here,  with  the  primary  consideration  being  given 
to  delta  vee.  The  time-of-fl ight  could  just  as  easily  have  been  the 
limiting  constraint  and  would  be  a  very  important  consideration  for  a 
manned  flight.  A  manned  flight  would  naturally  include  a  return  tra¬ 
jectory  and  consideration  would  be  required  for  how  long  to  remain  on 
the  asteroid.  A  trade-off  between  delta  vee  on  both  legs,  launch 
windows,  time-of-fl ight  and  time  on  the  asteroid  would  be  a  mission 
planning  task.  Mission  planning  is  beyond  the  scope  of  this  analysis. 
However,  the  tools  presented  in  these  computer  programs  provide  the  data 
that  would  be  required  for  this  kind  of  mission  planning. 

TABLE  II 

Favorable  Launch  Opportunities 
to  Asteroid  1982DB  from  1985  to  2020 


Total 

True  Anomaly 
of  Earth 

True  Anomaly 
of  Asteroid 

Time- 

Delta  V 

at  Launch 

at  Arrival 

of-Fl ight 

Date 

( km/s) 

(deg) 

(deg) 

(days) 

3 

Feb 

1991  4.5 

2 

13 

Feb 

1991  4.6 

3 

17 

Dec 

2001  5.6 

33 

26 

Jan 

2002  4.5 

1 

5 

Feb 

2002  4.5 

2 

15 

Feb 

2002  4.6 

2 

7 

Mar 

2002  5 .6 

4i 

16 

Mar 

2004  5.8 

6' 

29 

Jan 

2011  4.5 

2' 

29 

Jan 

2011  4.5 

2' 

28 

Feb 

2011  5.6 

5i 

TABLE  III 


Favorable  Launch  Opportunities 
to  Asteroid  1982XB  from  1985  to  2020 


Date 


Total 
Delta  V 
(km/s) 


True  Anomaly  True  Anomaly 
of  Earth  of  Asteroid 

at  Launch  at  Arrival 

(deg)  (deg) 


Time- 
of-Fl ight 
(days) 


21  Nov 

1987 

7.1 

310 

190 

1  Dec 

1987 

6.0 

320 

80 

21  Dec 

1987 

5.3 

340 

80 

30  Jan 

1988 

7.0 

20 

160 

24  Dec 

1992 

5.3 

340 

100 

3  Jan 

1993 

5.4 

350 

150 

2  Feb 

1993 

7.0 

20 

160 

8  Dec 

1997 

6.0 

310 

70 

17  Jan 

1998 

5.7 

360 

140 

6  Feb 

1998 

7.0 

20 

160 

22  Nov 

2002 

7.1 

300 

190 

2  Dec 

2002 

6.0 

320 

80 

31  Jan 

2003 

7.0 

20 

160 

26  Nov 

2007 

7.1 

310 

180 

4  Feb 

2008 

7.0 

20 

160 

23  Nov 

2017 

7.1 

310 

190 

TABLE  IV 

Favorable  Launch  Opportunities 
to  Asteroid  Anteros  from  1985  to  2020 


True  Anomaly  True  Anomaly 
Total  of  Earth  of  Asteroid  Time- 
Delta  V  at  Launch  at  Arrival  of-Flight 
(km/s)  (deg)  (deg)  (days) 


Date 


True  Anomaly 

True  Anomaly 

Total 

of  Earth 

of  Asteroid 

Time- 

Delta  V 

at  Launch 

at  Arrival 

of-Fl ight 

Date 

( km/  s ) 

(deg) 

(deg) 

(days) 

TABLE  VII 


Favorable  Launch  Opportunities 
to  Asteroid  1985TB  from  1985  to  2020 


1 

Dec 

1935  1 

20 

Jan 

1986  1 

21 

Jan 

2003  1 

20 

Jan 

2005  1 

19 

Jan 

2009  1 

4 

Apr 

2010  1 

18 

Jan 

2013  1 

17 

Jan 

2017  1 

L3.3 

320 

111 

[4.9 

10 

121 

L4 .9 

10 

121 

L6.9 

10 

291 

L6.8 

10 

281 

[6.9 

80 

341 

[6.5 

0 

281 

L6.8 

0 

271 

VII.  Suggestions  and  Recommendations 


The  results  presented  in  the  previous  chapter  are  essentially  a 
demonstration  of  the  potential  of  the  programs  developed  to  analyze  the 
accessibility  of  Earth-approaching  asteroids.  The  analysis  process  is 
complete  in  that  it  proceeds  from  fundamental  information  about  the 
orbits  of  the  asteroids  under  consideration  through  the  various  stages 
to  the  launch  dates  for  possible  rendezvous  missions.  This  series  of 
computer  programs  is  a  first  attempt  to  model  the  analysis  process  from 
start  to  finish.  Many  refinements  of  the  programs  are  readily  apparent 
to  make  these  tools  truly  useful.  Some  suggestions  for  these  improvements 
are  presented  in  this  chapter. 

Orbit  Accuracy 

All  the  orbital  motion  calculated  in  these  programs  has  been  ideal 
Keplerian  two-body  motion.  In  the  short  term,  for  one  or  two  years,  this 
assumption  is  reasonably  valid,  especially  considering  that  the  initial 
orbit  of  newly-discovered  asteroids  is  usually  not  known  very  well  any 
way.  For  asteroids  that  have  been  observed  on  several  orbits,  and  for 
which  the  orbital  elements  have  been  very  accurately  determined, 
increased  computational  accuracy  of  the  orbital  motion  may  be  desiraDle. 
This  is  especially  important  for  planning  rendezvous  missions  to  be  con¬ 
ducted  tens  of  years  in  the  future.  Perturbations  in  the  orbits  of  the 
bodies  involved  can  be  significant,  and  of  course  grow  with  time.  Some 
orbits  are  perturbed  very  little  while  others  could  be  greatly  perturbed, 
depending  on  their  proximity  to  the  orbits  of  the  larger  planets.  The 


asteroids  that  most  closely  approach  the  Earth  are  going  to  be  the  most 
perturbed  by  the  Earth. 

Considerable  effort  will  be  required  to  modify  these  programs  to 
account  for  perturbations.  However,  since  most  of  the  asteroids  in 
question  never  venture  too  far  from  the  Earth,  a  very  reasonable  approxi¬ 
mation  could  be  obtained  by  considering  only  the  perturbing  effects  of 
Jupiter,  Mars,  Venus  and  the  Earth.  But  even  this  simplification  still 
leaves  a  challenging  six-body  orbital  mechanics  problem  which  would 
probably  have  to  be  solved  by  numerical  methods  of  computer  iterations. 

This  is  a  slow  process  if  any  degree  of  accuracy  is  required.  Fortunately, 
the  current  computer  program  is  not  at  all  time  constrained  in  its  cal¬ 
culations  of  orbits.  Almost  all  of  the  15.5  minutes  of  computer  time 
required  by  the  main  program  is  tied  up  in  the  delta  vee  algorithms. 

The  orbits  are  calculated  completely  outside  of  this  time  consuming  part 
of  the  program.  There  would  be  no  multiplicative  effect  on  time  by 
enhancing  the  program  to  account  for  perturbations. 

Delta  Vee  Algorithms 

Most  of  the  time  required  to  solve  the  delta  vee  algorithms  is  con¬ 
sumed  in  an  inefficient  routine  that  finds  the  minimum  value  on  the 
curve  defined  by  the  two  equations  that  use  the  parameter  of  the  inter¬ 
cept  ellipse  as  a  free  variable.  About  30  iterations  of  a  bisection 
search  routine  are  required  to  converge  to  a  reasonable  accuracy.  The 
procedure  called  Kepler  in  the  program  CalcAnomal ies  was  tried  as  a 
replacement  for  the  bisection  routine.  Although  the  Kepler  procedure  is 
extremely  fast,  it  did  not  often  find  the  global  minimum  on  the  curve  but 
converged  on  a  local  minimum  and  was  abandoned.  However,  it  is  virtually 


guaranteed  to  find  a  solution  in  three  iterations.  If  the  problem  of 
local  convergence  could  be  overcome,  the  time  required  for  the  whole 
program  could  be  reduced  from  15.5  minutes  down  to  about  1.5  minutes. 
This  would  offset  the  inevitable  increase  in  time  that  would  result  if 
the  program  was  modified  to  account  for  perturbations. 


Integrating  All  the  Computer  Programs 


Except  for  the  validation  program  called  Hel ioCentricOrbit  which  is 
not  essential  to  the  analyses,  all  the  other  programs  could  easily  be 
combined.  They  were  left  separated  merely  as  a  convenience  during 
program  development  and  testing.  A  master  program,  using  menu-driven 
options,  would  relieve  most  of  the  manual  entries  required  by  the  present 
programs.  All  programs  rely  on  a  common  set  of  data  files  containing 
the  orbital  elements.  There  would  be  less  need  to  generate  inter¬ 
mediate  results  in  text  files  as  is  now  required.  The  process  of  com¬ 
bining  these  programs  has  been  facilitated  by  ensuring  that  all  of  the 
programming  is  modular  and  most  of  it  is  portable  to  other  programs. 


Launcn  Date  Search  Routi nes 

The  program  Fi ndlaunchDates  does  a  very  efficient  search  routine 
for  selected  combinations  of  true  anomalies  and  time-of-fl ight.  However, 
all  the  input  must  now  be  manually  entered  from  the  tabulated  values  of 
the  Prime  Rib  plots  and  time-of-fl ight  tables.  This  is  a  reasonable 
approach  if  only  the  launch  opportunities  for  a  limited  number  of  com¬ 
binations  of  true  anomalies  are  required.  Unfortunately,  this  does 
not  ensure  that  the  optimum  or  even  the  better  launch  opportunities  for 
a  given  period  are  found.  There  may  indeed  be  no  match  found  to  satisfy 


the  conditions  of  the  global  optimal  point  on  the  Prime  Rib  curves.  It 
then  becomes  a  matter  of  good  luck  or  laborious  manual  input  to  search 
for  the  next  best  solution.  There  are  1296  possible  inputs  making  an 
automated  global  search  a  very  desirable  feature.  Such  a  search  could 
be  conveniently  added  when  the  separate  programs  are  combined  into  a 
master  program. 

Detailed  Trajectories 

When  the  optimal  delta  vees  are  found  for  each  combination  of  true 
anomalies,  very  little  is  actually  known  about  the  actual  trajectory. 
However,  once  the  parameter  and  the  time-of-fl ight  of  the  trajectory  are 
calculated,  all  the  other  elements  of  the  ellipse  fall  out  of  the  equa¬ 
tions.  A  pictorial  or  tabular  representation  of  the  rendezvous  trajectory 
could  then  be  constructed  to  better  visualize  and  understand  the  geometry 
of  the  optimal  solutions.  This  information  could  reveal  fly-by  oppor¬ 
tunities  for  other  asteroids  enroute  to  the  main  objective.  It  could 
also  be  used  for  a  checx  on  potential  perturbations  by  observing  the 
path  of  the  trajectory  compared  to  the  positions  of  the  major  planets. 

Return  Mi ssions 

The  current  system  of  computer  programs  permit  return  missions  to 
be  analyzed  but  the  process  requires  a  considerable  amount  of  manual 
intervention.  The  addition  of  a  return  mission  search  routine  would 
greatly  enhance  the  utility  of  the  programs.  A  further  data  file  would 
be  required  containing  the  gravitational  constants  of  the  planets  and 
the  larger  asteroids.  The  current  programs  have  this  information  for 
only  the  Earth  and  a  representative  value  for  the  small  asteroids 


assuming  a  body  of  one  kilometer  in  diameter.  To  save  on  the  delta  vee 
requirements  for  return  missions,  some  are  proposed  which  use  atmos¬ 
pheric  braking  to  return  to  a  Shuttle-accessible  low  Earth  orbit.  A 
capability  to  handle  this  type  of  trajectory  would  be  useful  for  plan¬ 
ning  return  trajectories. 

Adjustable  Prime  Rib  Plots 

The  current  procedure  that  plots  the  Prime  Ribs  uses  a  single  fixed 
grey  scale  to  represent  the  delta  vee  values.  This  is  quite  useful  for 
comparing  several  plots  but  within  one  plot,  it  would  be  better  to  be 
able  to  adjust  the  grey  scale.  For  example,  it  might  be  desirable  to 
always  have  the  optimal  area  the  darkest  black,  regardless  of  the  actual 
delta  vee  that  it  represents.  It  might  be  practical  to  display  only  that 
part  of  the  plot  with  values  of  delta  vee  up  to  a  certain  value,  such  as 
the  maximum  delta  vee  available  in  a  proposed  rocket  system  for  a  par¬ 
ticular  mission. 

Conclusion 

The  series  of  programs  that  has  been  presented  in  this  thesis  pro¬ 
vides  a  basis  of  meaningful  analyses  of  the  accessibility  of  Earth- 
approaching  asteroids.  These  analyses  have  been  conducted  for  selected 
asteroids  and  the  results  are  apparent.  However,  the  potential  of  these 
programs  is  to  provide  a  more  accurate  and  more  usable  system  that  can 
be  implemented  on  modest  computer  resources  to  extend  the  capability  of, 
and  access  to,  these  kinds  of  analyses. 


Appendix  A 


PASCAL  Program  to  Implement  the  Hulk  ewe r- La u -Bender  A1 gorithms 


program  TrueAnomDeltaVee; 

(*  Calc  Delta  V  in  Space  of  True  Anom  +  Plot.  Obj  pos'n,  velocity  in  *) 
(*  heliocentric/ecliptic  I JK  ref  frame.  Then  runs  DeltaVee  program  on  *) 
(*  360x360  combinations  of  Departure/Arrival  True  Anom.  Hulkower,  Lau,*) 
(*  Bender  Method.  Plots  Prime  Ribs,  generates  data  files  for  Delta  Vee*) 
(*  and  Times-of-Fl i ght.  *) 


const 

MaxRecSize  =  25; 

Hel i oGravConst  =  3 . 9640157489E- 14 ;  (*  AU3/s2  derived  from 


Astronomical  Almanac  85  p. 

K6  *) 

AU  =  1.49597870E8; 

(*  km  deri ved  from 

Astronomical  Almanac  85  p. 

K6  *) 

GravConstSun  =  132718E6; 

(*  km3/s2  *) 

GravConstEarth  =  398603.0; 

(*  km3/s2  * ) 

type 

Name  =  string  [10]; 
ElementSet  =  record 

Asteroi dName: 

Name; 

(*  object  name  *) 

Tp: 

real ; 

(*  date  of  perih  passage  or  data 

pt  *) 

i  : 

real ; 

(*  inclination  to  ecliptic  *) 

w  : 

real ; 

(*  arguement  of  perihelion  *) 

N  : 

real ; 

(*  longitude  of  ascending  node  *) 

a  : 

rea  1 ; 

(*  semi -major  axis,  in  AU  *) 

e  : 

rea  1  ; 

(*  eccentri  city  *) 

Mo: 

real ; 

(*  mean  anomoly  at  To, 

Mo  =  0  at  perihel i on 

*) 

var 

AsteroidFile:  text; 

FileName:  string[14]; 

ElementSetRec:  array  [1. .MaxRecSize]  of  ElementSet; 
Ifor:  integer; 

MaxNoOfAsteroids:  integer; 

Selected:  integer;  (*  record  selected  to  work  with  *) 


Asteroi dName:  name; 

Tp:  real ; 
i  :  real; 
w  :  real ; 
N  :  real; 
a  :  real; 
e  :  real ; 


(*  date  of  perih  passage  or  data  pt  *) 
(*  inclination  to  ecliptic  *) 

(*  arguement  of  perihelion  *) 

(*  longitude  of  ascending  node  *) 

(*  semi-major  axis,  in  AU  *) 

(*  eccentri city  *) 


Mo:  real ; 


(*  mean  anomoly  at  To, 

Mo  =  0  at  perihelion  *) 

TrueAnomaly,  Eccentri cAnomaly  :  real; 

Isun,  Jsun,  Ksun:  real;  (*  heliocentric/ecliptic  coordinates  *) 
VelocI,  VelocJ,  VelocK:  real;  {*  heliocentric/ecliptic  velocity  *) 
Rsun:  real;  (*  distance  from  sun  *) 

A1,A2,B1  ,B2,Y1,Y2:  real;  (*  auxi 1 iaries  *) 

PiSun,  PjSun,  PkSun,  QiSun,  QjSun,  QkSun:  real;  (*  auxi 1 iaries  *) 
sini,  sinw,  sinN:  real;  (*  to  precompute  sin(i)  etc.  *) 
cosi ,  cosw,  cosN:  real; 

DeparturePositi on:  array [0.. 36,1. .3]  of  real; 

DepartureVelocity :  array  [0. .36,1 . .3]  of  real; 

Arri valPositi on:  array [0. .36,1 . .3]  of  real; 

Arri val Velocity :  array [0. .36,1 . .3]  of  real; 

Ul:  array[1..3]  of  real; 

U2:  array [1.. 3]  of  real; 

rl,r2:real ; 

DepartVel:  array [1.. 3]  of  real;  (*  old  coord  sys  *) 

ArrivalVel:  array[1..3]  of  real;  (*  old  coord  sys  *) 

DepartPos:  array [1.. 3]  of  real;  (*  old  coord  sys  *) 

ArrivalPos:  array [1. .3]  of  real;  (*  old  coord  sys  *) 

RlxR2:  array[1..3]  of  real; 

MagRlxR2  :  real; 

AngleBetweenVectors:  real; 
q  :  integer; 
counter:  integer; 

TrueEarth,  TrueTarget  :  integer; 

DataOut  :  text; 

OutFile  :  string[8]; 

DeltaVeeStore  :  array [0. .36,0. .36]  of  real; 

TimeOfFltStore  :  array [0. .36,0. .36]  of  real; 


function  power'a ,b: real ): real ; 
begi  n 

if  (a  =  0)  and  (b  =  0)  then  pcwer  :=  1 
else  begin 

if  a>0  then  power  :=  exp(b*ln(a)) 
else  begin 

if  a<0  then  begin 

writeln('  illegal  arg.  minus  power  '); 
power  :=  -999999999.9; 
end 


end; 


else  power  :=  0; 


end; 


end ; 


function  tan(x:  real)  :  real; 
begin 

tan  :=  sin(x)/cos(x); 
end;  (*  tan  *) 

function  arccos (a : real ) :  real; 
var 

temp:  real; 
begi  n 

if  a  >  1.0  then  a  :=  1.0; 
if  a  <  -1.0  then  a  :=  -1.0; 

if  a  =  0.0  then  a  :=  0.0000000001 ;  (*  avoid  divide  by  zero  *) 

temp  :=  arctan(sq  rt(1.0  -  a  *  a)  /  a); 

if  temp  <  0.0  then  arccos  :=  temp  +  Pi  else  arccos  :=  temp; 
end; 


procedure  GetDataFile;  (*  read  data  into  global  variables  *) 
begin 

Assign  (AsteroidFile, Filename);  (*  FileName  from  proc.  Init  Display  *) 
Reset  (AsteroidFile); 

Ifor  :=  0; 

while  not  eof  (AsteroidFile)  do 
begi  n 

ifor  :=  Ifor  +  1; 
wioh  ElementSetRecCl f or]  do 
begi  n 

readln  (AsteroidFile, 

AstercidName, 

TP, 

'  » 
w  , 

N  , 
a  , 
e  , 

Mo); 

end;  (*  with  *) 
end;  (*  while  *) 

MaxNoOf Asteroids  :=  Ifor; 

writeln( 'Objects  in  file'); 
writeln( 

'No.  Object  Tp  i  w  N  a  e  Mo' 


for  Ifor  :=  1  to  MaxNoOfAsteroids  do 
begin 

with  ElementSetRec[Ifor]  do 

writoln(AsteroidName: 10,Tp:13:4,i:9:4,w:9:4,N:9:4,a:8:4,e:7:4, 

Mo:9:4) ; 

end;  (*  for  *) 
writeln; 

end;  (*  GetOataFile  *) 


procedure  SelectObject ;  (*  Global  variables  *) 
begi  n 

write(  'Select  one  ...  '); 
readl n (Selected ) ; 
writeln; 

AsteroidName  :=  Element.SetRec[Selected]. AsteroidName; 

Tp  :=  E1ementSetRec[Selected].Tp; 
i  :=  ElementSetRec[Selected].i  ; 
w  :=  ElementSet.Rec[Selected].w  ; 

N  :=  ElementSetRec[Selected].N  ; 
a  :=  ElementSetRec[Selected].a  ; 
e  :=  ElementSetRec[Selected].e  ; 

Mo  :=  E1ementSetRec[Selected].Mo; 

writeln (AsteroidName.Tp: 10: 1 ,i :8:3,w:8:3,N:8:3,a:8:3,e: 8: 3, Mo: 8: 3) 
delay (20) ; 

end;  (*  SelectObject  *) 


procedure  InitialDisplay ; 
type 

FileNames  =  string[14]; 
var 

DataFi leSelecti on:  integer; 
FilesOnDisk:  array[1..6]  of  FileNames; 

begin 

F i lesOnDisk [1]  :=  'PLANET. ELE ' ; 
FilesOnDisk [2]  :=  ‘COMET. ELE'; 

F i lesOnDisk [3]  :=  ' AST-0001. ELE ' ; 
FilesOnDisk [4]  :=  'AST-0021. ELE ' ; 
FilesOnDisk [5]  :=  'AST-0041. ELE 1 ; 
FilesOnDisk [6]  :=  'ANTEROS.ELE ' ; 

cl  rscr ; 
writeln; 
writeln; 
writeln: 


writelnC  HELIOCENTRIC/ELLIPTIC  COORDINATES 

writelnC  for  360  deg.  of  TRUE  ANOMALIES  '); 

wri  teln; 

*riteln(‘  planets'); 

writelnC  comets'); 

writelnC  asteroids'); 

writeln; 

writelnC  Data  files:  contain  classic  orbit  elements.'); 
writelnC  Input:  Increment  for  TRUE  ANOMALY.'); 

writelnC  Output:  Hel i ocentri c/Ecl ipti c  Coordinates. ') ; 

writeln; 

writelnC  Data  Files  on  Disk . '); 

writelnC'  1.  PLANET. ELE ') ; 

writeln('  2.  COMET. ELE1); 

writelnC  3.  AST-0001. ELE ') ; 

writelnC  4.  AST-0021. ELE ') ; 

writelnC  5.  AST-0041. ELE ') ; 

writelnC  6.  ANTEROS.ELE  ' ) ; 

writeln('  7.  User-supplied'); 

writeln; 

write(‘  Select  (  1..7  )  '); 


DataFi leSelection  :=  0; 
while  not  (DataFi leSelecti on  in  [1..7])  do 
readln (DataFi leSelection); 
if  DataFi leSelecti on  <  7  then  Filename  := 

FilesOnDisk  [DataFileSelection] 

else 

begin 

writeln; 

write( 'Input  FILENAME. EXTENTION  '); 
readln (Filename) ; 
writeln; 
end; 

end;  (*  Initial  Display  *) 

procedure  ConvertToRadians ; 
begin 

w  :=  w  *  Pi  /  180.0; 

i  :=  i  *  Pi  /  180.0; 

N  :=  N  *  Pi  /  180.0; 

end; 

procedure  PreComputeT ranscendental s ; 
begi  n 

sini  :=  sin(i);  cosi  :=  cos(i); 
sinw  :=  sin(w);  cosw  :=  cos(w); 
sinN  :=  sin(N);  cosN  :=  cos(N); 


procedure  Cal cAuxQuant  ; 
begin 


A1  :=  sinN*sinw;  (*  auxiliary  quantities  *) 

B1  :=  cosN*sinw; 

Y1  :=  sini*sinw; 

A2  :=  sinN*cosw; 

B2  :=  cosN*cosw; 

Y2  :=  sini*cosw; 

(*  GAUSSIAN  CONSTANTS  *) 

PiSun  :=  B2  -  Al*cosi ; 

PjSun  :=  (A2  +  Bl*cosi); 

Pk Sun  :=  Yl; 

QiSun  :=  -B1  -  A2*cosi; 

QjSun  :=  (-A1  +  B2*cosi); 

Qk  Sun  :=  Y2; 

end;  (*  CalcAuxQuant  *) 


procedure  Cal culateObjPositi on; 

(*  Isun,  etc  are  in  AU's  *) 
begi  n 

Eccentri cAnomaly  :=  2.0  *  arctan(  tan(0.5  *  TrueAnomaly  *  Pi/130.0) 

*  sq  rt  ( ( 1 . 0  -  e)  /(1.0  +  e))  ); 

Isun  :=  a *Pi Sun*(cos (Eccentri cAnomaly  )-e)  + 

a*sq  rt (1.0-e*e )*Qi Sun *s in (Eccentri  cAnomaly  ) ; 

Jsun  :=  a*PjSun*(cos (Eccentri cAnomaly  )-e )+ 

a*sq  rt (1.0-e*e )*Qj Sun *s in (Eccentri cAnomaly  ) ; 

Ksun  :=  a*Pk Sun*(cos (Eccentri cAnomaly  )-e )+ 

a*sq  rt(1.0-e*e)*0kSun*sin(EccentricAnomaly  ); 

Rsun  :=  sq  rt (Isun*Isun  +  Jsun*Jsun  +  K.sun*Ksun);  (*  in  AU  *) 

end;  (*  Cal culateObjPositi on  *) 


procedure  Cal culateObjVel ocity ; 

(*  Velocities  are  in  AU's  *) 
var 

parameter:  real; 

sq  rtHp,sinT,cosT:  real;  (*  temporaries  *) 


IV  V 


sinT  :=  sin (TrueAnomaly  *  Pi/180.0); 
cost  :=  cos (TrueAnomaly  *  Pi/180.0); 

parameter  :=  a  *  (l.C  -  e  *  e); 

sqrtHp  :=  sq  rt (Hel i oGravConst  /  parameter)  ; 


VelocI  :=  sqrtHp  *  (  (-sinT  *  PiSun)  +  (e  +  cosT)  *  QiSun  ) 

VelocJ  :=  sqrtHp  *  (  (-sinT  *  PjSun)  +  (e  +  cosT)  *  QjSun  ) 

VelocK  :=  sqrtHp  *  (  (-sinT  *  PkSun)  +  (e  +  cosT)  *  QkSun  ) 

end;  (*  CalculateObjVelocity  *) 


procedure  Departure360; 
var 

counter:  integer; 

begi  n 
writeln; 
writeln; 

CalcAuxQuant; 

TrueAnomaly  :=  0.0; 

while  TrueAnomaly  <=  360.0  do 
begi  n 

counter  :=  round(TrueAnomaly /10.0) ; 

Cal culateObjPositi on; 

CalculateObjVelocity; 

DeparturePositi on  [counter, 1]  :=  Isun  *  All 
DeparturePositi on  [counter, 2]  :=  Jsun  *  All 
DeparturePosition  [counter, 3]  :=  Ksun  *  AU 


DepartureVel oci ty  [counter, 1]  :=  VelocI  *  AU; 
DepartureVel oci ty  [counter, 2]  :=  VelocJ  *  AU; 
DepartureVelocity  [counter, 3]  :=  VelocK  *  AU; 

TrueAnomaly  :=  TrueAnomaly  +  10.0; 
end;  (*  while  *) 
end;  (*  Departure360  *) 

procedure  Arrival 360; 
var 

counter:  integer; 


begin 

writeln; 

writeln; 

Cal cAuxQuant; 

TrueAnomaly  :=  0.0; 

while  TrueAnomaly  <=  360.0  do 


A-7 


counter  :=  round (TrueAnomaly /10.0) ; 

Cal culateObjPositi on ; 

Cal  culateOb.jVelocity ; 

Arri valPosition  [counter, 1]  :=  Isun  *  AU  ; 
Arri valPositi on  [counter, 2]  :=  Jsun  *  AU  ; 
Arri valPosition  [counter, 3]  :=  Ksun  *  AU  ; 


Arri val Vel ocity  [counter.1]  :=  VelocI  *  AU 
Arrival  Velocity  [counter, 2]  :=  VelocJ  *  AU 
Arri val Velocity  [counter, 3]  :=  VelocK  *  AU 

TrueAnomaly  :=  TrueAnomaly  +  10.0; 
end ;  (*  while  *) 
end;  (*  Arri val 360  *) 


procedure  RlcrossR2; 
begi  n 

RlxR2['l]  :=  DepartPos[2]  *  Arri valPos[3] 

-  DepartPos[3]  *  Arri valPos[2]  ; 

RixR2[2]  :=  DepartPos[3]  *  Arri valPos[l] 

-  DepartPos[l]  *  Arri valPos[3]  ; 

RlxR2[3]  :=  DepartPos[l]  *  Arri valPos[2] 

-  DepartPos[2]  *  Arri valPos[l]  ; 


MagRlxR2  :=  sqrt(  RlxR2[l]  *  RlxR2[l] 

+  RlxP.2[2]  *  RlxR2[2] 

+  RlxR2[3]  *  RlxR2[3]  )  ; 

end;  (*  RlcrossR2  *) 


procedure  AngleRlR2: 
var 

hearth:  array[1..3]  of  real: 
temp:  real; 

begi  n 

AngleBetweenVectors  :=  arccos(  (  DepartPos[l]  *  Arri valPos[l] 

+  DepartPos[2]  *  Arri valPos[2] 

+  DepartPos[3]  *  Arri valPos[3]  ) 

/  (  rl  *  r2  )  ); 


RlcrossR2; 


(*  cal culate  q  *) 


hearth[l]  :=  DepartPos[2]  *  DepartVel[3] 

-  DepartPos[3]  *  DepartVel [2]  ; 

hearth[2]  :=  DepartPos[3]  *  DepartVel[l] 

-  Depart.Pos[l]  *  DepartVel [3]  ; 

hearth[3]  :=  DepartPos[l]  *  DepartVel [2] 

-  DepartPos[2]  *  DepartVel [1]  ; 


temp  :=  RlxR2[l]  *  hearth[l] 

+  RlxR2[2]  *  hearth[2] 

+  RlxR2[3]  *  hearth[3]  ; 

if  temp  >=  0.0  then  q  :=  1  else  q  :=  -1  ; 


If  q  =  -1  then  AngleBetweenVectors  :=  2.0  *  Pi 

-  AngleBetweenVectors; 

end;  (*  AngleRlR2  *) 


procedure  GetNextSetT rueAnomal ies ; 
var 

counter  :  integer; 
begin 

for  counter  :  =  1  to  3  do 
begin 

DepartPos[counter]  :=  DeparturePcsiti on[TrueEarth, counter]; 
Arri va!Pos[counter]  :=  Arri va!Positior[TrueTarget, counter]; 

DepartVel [counter]  :=  DepartureVel ocity  [TrueEarth, counter]; 

Arri val Vel [counter]  :=  Arri val Vel ocity [TrueTarget, counter]; 
end; 

ri  :=  sq  rt (DepartPos[l]  *  DepertPos[l] 
hDepartPos[2]  *  DepartPos[2] 
fDepartPos[3]  *  DepartPos[3]) ; 

r2  :=  sq  rt(Arri valPos[l]  *  Arri val Pos[l] 

+Arri valPos[2]  *  Arri valPos[2] 

+Arri valPos[3]  *  Arri val Pos[3]) ; 

U 1  [ 1 ]  :=  OeparturePosition  [TrueEarth, 1]  /  rl; 

U 1  [2]  :=  OeparturePosition  [TrueEarth, 2]  /  ri; 

Ul[3]  :=  OeparturePosition  [TrueEarth,3]  /  ri; 

U2[l]  :=  Arri val Position  [TrueTarget,l]  /  r2; 

U2[2]  :=  Arri val Positi on  [TrueTarget ,2]  /  r2; 

U2[3]  :=  ArrivalPosition  [TrueTarget, 3]  /  r2; 

end;  (*  proc  GetNextSetTrueAnomalies  *) 
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procedure  Shade(X,Y:  integer;  c 
var 

s,a,b,d,e:  integer; 
begin 

if  c  <=  4.5  then 
begin 

for  s  :=  0  to  2  do 
begin 

for  a  :=  0  to  3  do 
begi  r, 

for  b  :=  0  to  3  do 
begi  n 

plot(x+s*4+a ,y+b,l) ; 
end; 

end; 

end; 

end; 


if  (c  <=  5.0)  and  (c  >  4.5)  then 
begi  n 

for  s  :=  0  to  2  do 
begin 

plot (x+s *4+1, y  ,1 ) ; 
p lot (x+s *4+2, y  ,  1 ) ; 
plot(x+s*4+3,y  ,1); 
pi ot (x+s*4,y +  1 , 1 ) ; 
p  1  ot ( x+s  *4+2  ,y  + 1 , 1 ) ; 
plot(x+s *4+3, y  + 1,1 ) ; 
plot. (x+s *4,y+2,l ) ; 
pi  ot  (x  +s*4+l  ,y +2 , 1 ) ; 
plot  (x+s*4+3,y+2,l ) ; 
plot(x+s*4,y+3,l); 
plot (x+s *4+l,y+3,l); 
pi  ot  (x+s *4+2  ,y  +3, 1 ) ; 
end ; 

end; 


if  (c  <=  6.0)  and  (c  >  5.0)  then 
begi  n 

for  s  :=  0  to  2  do 
begi  n 

pi ot (x+s *4 ,y  ,  1 ) ; 
plot(x+s*4+3,y  ,1 ) ; 
pi  ot  (x+s*4+l,y  +  l  ,1 ) ; 
plot(x+s *4+2,y +  1 , 1 ) ; 
plot(x+s*4+l,y+2,l); 
plot(x+s*4+2,y+2,l); 
pl  ot  (x+s*4,y +3,1 ) ; 


p1ot(x+s*4+3,y+3,l) ; 
end; 

end; 

if  (c  <=  8.0)  and  (c  >  6.0)  then 
begi  n 

for  s  :=  0  to  2  do 
begin 

plot(x+s*4,y  ,1); 
plot (x+s  *4+3  ,y  ,  1 ) ; 
plot (x+s  *4+l,y  +  l,l); 
p  1  ot  ( x+s  *4+2  ,y  +2 , 1 ) ; 
pi  ot  (x+s*4,y+3,l ) ; 
plot  (x+s *4+3 ,y+3,l ) ; 
end; 

end ; 

if  (c  <=  10,0)  and  (c  >  8.0)  then 
begi  n 

for  s  :=  0  to  2  do 
begi  n 

plot  (x+s  *4,y  ,1); 
plot(x-ks*4+2,y+l,l); 
plot  (x+s *4+l,y+2,l); 
plot (x+s *4+2, y +4,1 ) ; 
end; 

end; 

if  (c  <=  14.0)  and  (c  >  10.0)  then 
begin 

for  s  :=  0  to  2  do 
begin 

pi  ot  (x+s  *4+1  ,y-rl  ,1 ) ; 
pi ot (x+s *4+2, y +2,1 ) ; 
end ; 
end ; 

if  (c  <  20.0)  and  (c  >  14.0)  then 
begin 

for  s  :=  0  to  2  do 
begi  n 

plot (x+s *4+2, y +2,1 ) ; 
end ; 

end ; 

end;  (*  Shade  *' 

procedure  Cal  cl ; 
var 

V  :  array [1. .3]  of  real ; 
vl  :  array[1..3]  of  real; 
v2  :  array [1.. 3]  of  real; 
pi  :  array [1. .3]  of  real ; 


! 


p2  :  array [1.. 3]  of  real; 

1 1  :  array  [1. .3]  of  real  ; 

12  :  array [1. .3]  of  real ; 
dll  :  array[1..3]  of  real; 
dl 2  :  array[1..3]  of  real; 


dVl,dV2  :  real 
z  :  real 
p,pn,pnl  :  real 
pmin,  pmax  :  real 
sq  rtupt  :  real 
SL  :  real 
TempQql,  TempQq2  :  real 
rldotr2  :  real 
Ildotdll  :  real 
I2dotdI2  :  real 
TempNumer,  TempDenom  :  real 
Theta  :  real 
Theta360  :  real 
t  :  real 
J,  TernpJ  :  real 
sq  rtup  :  real 
Ildotll  :  real 
I2dotI2  :  real 
Diff  :  real 
Tempp  :  real 


•  «  ’•  rm  -  v 


const 

u  =  132718E6 ;  (*  km3/s2  *) 

q 1  -  6656.3; 
q  2  =  15.0; 

QQ1  =  6656.0; 

QQ2  =  5E1 ; 

uul  =  3.936005E5;  (*  km3/s2  *) 

uu2  =  3.37E-7;  (*  kn3/s2  Direct  proportion  to  volume  1982DB 

/  volume  Earth  *) 
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Procedure  FindDiff; 
begi  n 

sq  rtup  :=  sq  rt(u  *  p ) ; 

z  :=  sqrt(u  /  p)  *  tan(Theta/2.0) ; 
for  counter  :=  1  to  3  do 
begin 

v[counter]  :=  sq rtup  *  (p2[counter]  -  pl[counter]  )  /  MagRlxR2; 


Il[counter]  := 
I2[counter]  := 
dIl[counter]  := 
dI2[counter]  := 
end; 


-  vl[counter]  +  SL  *  (  v[counter]  +  Z  *  Ul[counter]); 
v2[cou nter]  -  SL  *  (  v[counter]  -  Z  *  U2[counter]) ; 
SL  *  0.5  /  p  *  (v[counter]  -  l  *  Ul[counter]) ; 

-SL  *  0.5  /  p  *  (v[counter]  +  z  *  U2[counter]) ; 
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,  '*  V,4* fcV  »'* 


TempQql  :=  sqrt(2.0  *  uul  *  QQ 1  / (q  l*(QQl-tq  1 ) )  )  ; 

TempQq  2  :=  sqrt(2.0  *  uu2  *  QQ2/(q  2*(QQ2*q  2) )  )  ; 

IldotU  :=  T  1[1]  *  Il[l]  +  1 1  [ 2 ]  *  1 1[2]  +  1 1  [3 ]  *  Il[3]; 

1 2dot 1 2  :=  I2[l]  *  I2[l]  +  I2[2]  *  I2[2]  +  I2[3]  *  I2[3] ; 

dVl  :=  sqrt(  Ildotll  +  (2.0  *  uul/ql)  )  -  TempQql  ; 

dV2  :=  sqrt(  I2dotI2  +  (2.0  *  uu2/^  2 )  )  -  TempQq 2  ; 

Ildotdll  :=  Il[l]  *  d 1 1 [ 1 ]  +  1 1[2]  *  dll[2]  +  Il[3]  *  dll[3]; 
I2dotdI2  :=  I2[l]  *  dI2[l]  +  I2[2]  *  dI2[2]  +  I2[3]  *  dl 2[3] ; 

Oiff  :=  Ildotdll  /  (  dVl  +  TempQql  )  *  I2dotdI2 

/  (  dV2  +  TempQq2  )  ; 

end;  (*  FindDiff  *) 


procedure  Newton; 
var 

dVee,  dVeen,  k  :  real ; 

Deltap  :  real ; 
dVeeOld  :  real; 

c  :  integer; 

const 

accuracy  =  0.0001; 
begi  n 

Deltap  :=  10.0; 

Dvee  :=  100.0; 
dVeeOld:-  1000.0; 

p  :=  (pflax  +  pMin)/2.0; 

k  :=  6.5  *  p; 
c  :=  0; 

while  (abs(dVee-dVeeOld)  >  accuracy)  do 
begi  n 

c  :  =  c  +  1 ; 
dVeeOld  :=  dVee; 

FindDiff; 
dVee  :=  dVl+dV2; 
p  :=  p+Deltap; 

FindDiff; 

dVeen  :=  dVl+dV2; 

if  dVeen  >  dVee  then 


if  c  >  50  then  dVee  :=  dVeeOld; 
end; 

J  :=  d V 1  +d V 2 ; 
end;  (*  Newton  *) 


procedure  GaussGodalTime0fFlight(rl,r2, Theta, p:  real;  var  t:  real); 
const 


u  = 

123718E15; 

(*  m3/s2 

var 

sinHalfTheta  : 

real ; 

cosHalfTheta  : 

rea  1 ; 

cosz 

.sinz  : 

real ; 

z 

real ; 

A,  B 

real ; 

Bmcosz  : 

real ; 

begin 

rl  :=  rl 

*  1000.0; 

r2  :=  r2 

*  1000.0; 

p  :=  p 

*  1000.0; 

SinHalfTheta  :=  sin(Tneta/2.0)  ; 
cosHalfTheta  :=  cos (Theta /2.0)  ; 

cosz  :=  (-2.0  *  rl  *  r2  *  sinHalfTheta  *  sinHalfTheta/p  +  rl  +  r2) 

/2.0/sq  rt(rl)/sq  rt(r2) /cosHalfTheta  ; 

Bmcosz  :=  sq  rt  (rl  )*sq  rt  (r2)  *  sinHalfTheta  *  sinHalfTheta  /  p  / 
cosHalfTheta  ; 

z  :=  arcccs (cosz ) ; 
sinz  :=  sin(z  ) ; 

A  :=  2.0  /  sqrt(u)  *  power(rl  ,0.75)*power(r2,0.75)  * 

pcwer(abs(cosHalfTheta ),1.5)  ; 

B  :=  (rl  +  r2)  /  2.0  /  sqrt(rl)  /  sqrt(r2)  /  cosHalfTheta  ; 

t  :=  A  *  sq  rt  (a  bs  (Bmcosz  ) )  *  (1.0  +  (Bmcosz)  *  (2.0*z  -  sin(2.0*z))/ 

2.0  /  s i nz  /  s i nz  /  sinz  ) ; 

t  :=  abs (t  /  3600.0  /  24.0)  ; 
end;  (*  GaussGodelTimeOfFl i ght  *) 


begin 

for  counter  :=  1  to  3  do 
begin 

pl[counter]  :=  DepartPos[counter]; 
p2[counter]  :=  Arri valPos[counter]; 


vl[counter]  :=  DepartVel [counter]; 
v2[counter]  :=  Arri val Vel [counter] ; 
end; 

Theta  :=  AngleBetweenVectors  ; 

Theta360  :=  Theta; 

if  theta  >  Pi  then  theta  :=  2.0  *  Pi  -  theta; 

rldotr2  :=  pl[l]  *  p2[l]  +  pl[2]  *  p2[2]  +  pl[3]  *  p2[3] 

TempNumer  :=  rl  *  r2  -  rldotr2; 

TempDenom  :=  sqrt(  2.0  *  (rl  *  r2  +  rldotr2)  )  ; 

prnin  :=  tempNumer  /  (rl  +  r2  +  TempDenom); 
pmax  :=  tempNumer  /  (rl  +  r2  -  TempDenom); 

if  Theta360  >  Pi  then  SL  :=  -1.0  else  SL  :=  1.0  ; 

Newton; 

GaussGodalTimeOfFl i ght (rl , r2,  Theta 360, p,t ) ; 

DeltaVeeSt.ore[TrueEarth,TrueTarget]  :=  J; 
TimeOfFltStore[TrueEarth,TrueTarget]  :=  t; 

(*  writeln(10  *  TrueEarth: 7 ,10  *  TrueTarget:7,J:10:4, '  ' , t : 6 : 1 

shade(85  +  TrueTarget*12,45  +  TrueEarth*4,J } ; 


end;  (*  Ca 1  cl  *) 


procedure  DeltaV360; 
var 

count  :  integer; 
counter:  integer; 
begin 

for  count  :=  0  to  36  do 
begi  n 

TrueEarth  :=  count  ; 
for  counter  :=  0  to  36  do 
begi  n 

TrueTarget  :=  counter; 
GetNextSetT  rueAnoma lies; 
AngleRlR2 ; 

Cal  cl ; 

end;  (*  counter  *) 
end;  (*  count  *) 


end;  (*  proc  DeltaV360  *) 


procedure  InitOutFi le; 
begi  n 

writeln;  write('What  is  the  name  (prefix)  of  the  output  file 
readl n (OutFi le ) ; 

assi gn(DataOut,  OutFi le+‘ .OUT 1 ) ; 
rewrite(DataOut); 
end;  (*  InitOutFile  *) 

procedure  WriteDVFile; 
var 

i , page, column ,rcw  :  integer; 
begi  n 

for  page  :=  0  to  2  do 
begi  n 

write(DataOut, 'Arr®0ep ' ); 
for  i  :=  0  to  11  do 
begi  n 

w rite (DataOut , (i +12*page  )*10: 5) ; 
end;  (*  for  i  *) 
write In (DataOut); 

for  rcw  :=  36  downto  0  do 
begin 

write(Data0ut,rcw*10:3, ‘  '); 

for  column  :=  0  +  12*page  to  11  +  12*page  do 
begin 

write (DataOut, DeltaVeeStoreCcolumn, row]: 5: 1 ); 
end;  (*  for  column  *) 
writeln(DataOut); 
end;  (*  for  row  *) 

writeln(DataOut);  writeln(DataOut);  writeln(DataOut); 
end;  (*  for  page  *) 

end;  (*  WriteDVOut  *) 


procedure  WriteTOFFile; 
var 

i , page, column, row  :  integer; 
begin 

for  page  :=  0  to  2  do 
begi  n 

write (DataOut, 'Arr®Dep 1 ) ; 
for  i  :=  0  to  11  do 
begin 

write (DataOut, (i +12*page  )*10: 5) ; 
end;  (*  for  i  *) 
writeln(DataOut); 

for  row  :=  36  downto  0  do 
begin 


write(Data0ut,row*10:3, 1  '); 

for  column  :=  0  +  12*page  to  11  +  12*page  do 
begin 

write(Data0ut,Time0fFltStore[column,rcw]:5:0); 
end;  (*  for  column  *) 
writeln(DataOut); 
end;  (*  for  row  *) 

writeln(DataOut);  writel n{DataOut ) ;  writeln(DataOut) 
end;  (*  for  page  *) 


end;  (*  WriteTOFOut  *) 


procedure  InitPrimeRibs; 


var 

i  :  integer; 
begi  n 

draw (  74,  193,  539,  193,  1) 
draw(  74,  44,  539,  44,  1) 

draw (  84,  197,  84,  40,  1) 


draw(529,  197,  529,  40,  1) 


for  i  :=  1  to  8  do 
begin 


draw(  84+48*i ,  193, 
draw(  84+48*i ,  44, 

84,  44+16*i , 

44+16*i , 
*) 


draw  ( 
draw(  529, 
end;  (*  for  i 


84+48*i , 
84+48*i , 

74,  44+16*i , 
539,  44+16*i , 


197, 

40, 


1) 

1) 

1) 

1) 


end;  (*  InitPrimeRibs  *) 


begin  (*  main  *) 

InitialDi splay ; 

GetDataFi le; 

SelectObject; 
ConvertToRadians ; 
PreComputeT  ranscendenta 1 s ; 
Departure360; 
InitialDisplay ; 

GetDataFi le; 

SelectObject; 
ConvertToRadians ; 
PreComputeT  rans  cendental s ; 
Arri val 360; 

InitOutFile; 

Cl rScr; 

HiRes; 

HiResCol or(7 ) ; 
InitPrimeRibs; 

DeltaV360; 


Appendix  B. 


PASCAL  Program  to  Implement 
the  Ross-Hulk  ower  A1 gorithms 


program  RossHulk owerOeltaVee; 


(*  Ross-Kulk ower  algorithm  to  find  optimal  delta  vee.  Calculates  *) 

(*  objects  position  and  velocity  in  heliocentric/ecliptic  UK.  ref  *) 
(*  frame.  Then  runs  OeltaVee  program  on  selected  combinations  of  *) 

(*  Departure/Arri  val  True  Anomalies.  *) 


const 

MaxRecSize  = 

Hel i oGravConst  =  3. 96401 57489E-14;  (*  AU3/s2  derived  from 

Astronomical  Almanac  85  p.  K6  *) 
AU  =  1.49597870E8;  (*  km  derived  from 

Astronomical  Almanac  85  p.  K6  *) 


type 

Name  =  string  [10]; 
ElementSet  =  record 
AsteroidName: 

Tp: 
i  : 
w  : 
N  : 
a  : 
e  : 
Mo: 


Name; 

(* 

object  name  *) 

real ; 

(* 

date  of  perihelion  passage/epoch  *) 

real ; 

(* 

inclination  to  ecliptic  *) 

real ; 

(* 

arguement  of  perihelion  *) 

real ; 

(* 

longitude  of  ascending  node  *) 

real  ; 

(* 

semi -major  axis,  in  AU  *) 

real ; 

(* 

eccentricity  *) 

real ; 

(* 

mean  anomoly  at  To, 

Mo  =  0  at  peri hel i on  *) 


end; 


var 

AsteroidFi le:  text; 

FileName:  string[14]; 

ElementSetRec:  array  [1. .MaxRecSize]  of  ElementSet; 
Ifor:  integer; 

MaxNoOfAsteroids:  integer; 

Selected:  integer;  (*  record  selected  to  work  with  *) 


AsteroidName:  name 
Tp:  real 
i  :  real 
w  :  real 
N  :  real 
a  :  real 
e  :  real 
Mo:  real 


(*  date  of  perihelion  passage/epoch  *) 
(*  inclination  to  ecliptic  *) 

(*  arguement  of  perihelion  *) 

(*  longitude  of  ascending  node  *) 

(*  semi -major  axis,  in  AU  *) 

(*  eccentri city  *) 

(*  mean  anomoly  at  To, 

Mo  =  0  at  perihelion  *) 


TrueAnomaly,  Eccentri  cAnomaly  :  real; 

Isun,  Jsun,  Ksun:  real;  (*  hel i ocentri c/ecl ipti c  coord  *) 


VelocI,  VelocJ,  VelocK:  real;(*  heliocentric/ecliptic  velocity  *) 
Rsun:  real;  (*  distance  from  sun  *) 

A1 ,A2 ,B1 ,B2 ,  Y1 , Y2:  real;  (*  auxi 1 iaries  *) 

PiSun,  PjSun,  PkSun,  QiSun,  QjSun,  QkSun:  real;  (*  auxiliaries  *) 
sini,  sinw,  sinN:  real;  (*  to  precompute  sin(i)  etc.  *) 
cosi ,  cosw,  cosN;  real; 

DeparturePosition:  array [0.. 35,1. .3]  of  real; 

DepartureVelocity :  array[0..35,1..3]  of  real; 

Arri valPositi on:  array [0. .35,1. .3]  of  real; 

Arri val Vel ocity  :  array [0. .35,1. .3]  of  real; 

i cPi rl,i cMi rl ,i cPi r2,i cMi r2,h:  array[1..3]  of  real; 

Uniti cPi rl,Uniti cMi rl.Uniti cPi r2,Uniti cMi r2:  array[1..3]  of  real; 

rl ,r2: real ; 

CC  :  real  ; 
q ,  SG  :  real  : 

PI, P2, Ml, M2:  real; 

DepartVelocOld:  array[1..3]  of  real;  (*  old  coord  sys  *) 

DepartVel ocNew:  array [1.. 3]  of  real;  (*  old  coord  sys  *) 

Arri veVel ocOl d:  array[1..3]  of  real;  (*  old  coord  sys  *) 

Arri veVel ocNew:  array[1..3]  of  real;  (*  old  coord  sys  *) 

DeparturePositOld:  array[1..3]  of  real;  (*  old  coord  sys  *) 

Arri val Posi tO I d:  array[1..3]  of  real;  (*  old  coord  sys  *) 

DeparturePositNew:  array [1.. 3]  of  real;  (*  new  coord  sys  *) 

Arri valPositNew:  array [1.. 3]  of  real;  (*  new  coord  sys  *) 

RlxR2:  array[1..3]  of  real; 

AngleBetweenVectors:  real; 

counter:  integer; 

TrueStar,  TrueP  :  integer; 


function  tan(x:  real)  :  real; 
begin 

tan  :=  sin(x)/cos(x); 
end;  (*  tan  *) 

(*  ***  The  following  8  procedures  are  identical  ***  *) 
(*  ***  to  procedures  of  the  same  names  in  the  *****  *) 
(*  ***  program  Hulk owerLauBenderDeltaVee  **********  *) 

procedure  GetDataFile;  (*  Reads  data  into  global  variables  *) 

procedure  SelectObject;  (*  Global  variables  *) 

procedure  InitialDisplay ; 
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procedure  ConvertToRadians ; 
procedure  PreC omputeT ranscendentals ; 
procedure  CalcAuxQuant  ; 
procedure  Cal culateObjPositi on; 
procedure  CalculateObjVel ocity ; 


procedure  Departure360; 


(*  This  procedure  is  modified  to  operate  on  a  single  true  *) 

(*  anomaly  rather  than  the  full  range  of  0-360  because  of  *) 

(*  the  instability  of  the  algorithms.  Problems  of  *) 

(*  convergence  may  occur  at  some  values  of  departure  true  *) 

(*  anomalies.  The  procedure  will  operate  through  the  0-360  *) 

i *  if  it  is  reconfigured  like  the  Departure360  algorithm  *) 

(*  in  the  program  Hulk owerLauBenderDeltaVee.  *) 


var 

counter:  integer; 

begin 

writeln; 

writeln; 

Cal cAuxQuant; 

TrueAricmaly  :=  110.0; 

while  TrueAnomaly  <=  110.0  do 
begi  n 

counter  :=  round(TrueAnomaly/10.0); 

Cal  culateObjPositi on; 

CalculateObjVelocity ; 

DeparturePositi on  [counter, 1]  :=  Isun; 

DeparturePositi on  [counter, 2]  :=  Jsun; 

DeparturePositi  on  [counter, 3]  :=  Ksun; 

DepartureVel oci ty  [counter, 1]  :=  VelocI  *  AU 

DepartureVelocity  [counter, 2]  :=  VelocJ  *  AU 

DepartureVel oci ty  [counter, 3]  :=  VelocK  *  AU 

TrueAnomaly  :=  TrueAnomaly  +  10.0; 
end;  (*  while  *) 
end;  (*  Departure360  *) 

procedure  Arrival 360; 
var 

counter:  integer; 
begin 
writeln; 
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writeln; 

Cal cAuxQuant ; 

TrueAnomaly  :=  0.0; 

while  TrueAnomaly  <=  360.0  do 

begi  n 

counter  :=  round (TrueAnomaly /10.0); 

Cal culateObjPositi on; 

Cal culateObjVel ocity ; 

Arri valPosition  [counter, 1]  :=  Isun; 

Arri valPositi on  [counter, 2]  :=  Jsun; 

Arri val Positi on  [counter, 3]  :=  Ksun; 

Arri val Vel ocity  [counter, 1]  :=  VelocI  *  AU 

Arri val Velocity  [counter, 2]  :=  VelocJ  *  AU 

Arri val Vel ocity  [counter, 3]  :=  VelocK  *  AU 

TrueAnomaly  :=  TrueAnomaly  +  10.0; 
end;  (*  while  *) 
end;  (*  Arri val 360  *) 


procedure  GetAl phaBetaGamma ; 
var 

hstar  :  array[1..3]  of  real; 

temp  :  real; 

begi  n 

DepartVel ocNew[l]  :=  OepartVel oc01d[l]  *  Uniti cPi rl[l] 

+  DepartVel oc01d[2]  *  Uniti cPi rl[2] 

+  DepartVel oc01d[3]  *  Uniti cPi rl[3] 

DepartVel ocNew[2]  :=  DepartVeloc01d[l]  *  Uniti cMiri[l] 

+  DepartVel oc01d[2]  *  Uniti cMi rl[2] 

+  DepartVeloc01d[3]  *  Uniti cMi rl[3] 

DepartVel ocNew[3]  :=  DepartVe!oc01d[l]  *  h [ 1 ] 

+  DepartVeloc01d[2]  *  h[2] 

+  DepartVelocQ?d[3]  *  h[3]  ; 


Arri veVel ocNew[l]  :=  Arri veVeloc01d[l]  *  Uniti cM i r 2 [ 1 3 

+  Arri veVel oc01d[2]  *  Uniti cMi r2[2] 
+  ArriveVeloc01d[3]  *  Uniti cMi r2[3] 

Arri veVel ocNew[2]  :=  Arri veVel ocOl d[l]  *  Uniti cPir2[l] 

+  Arri veVel oc01d[2]  *  Uniti cPi r2[2] 
+  Arri veVel oc01d[3]  *  Uniti cPi r2[3] 

Arri veVel ocNew[3]  :=  Arri veVeloc01d[l]  *  h[l] 

+  Arri veVel oc01d[2]  *  h[2] 

+  Arri veVel oc01ti[3]  *  h[3]  ; 


DepartVel ocNew[l]  :=  DepartVel ocNew[l]  *  PI 
OepartVel ocNew[2]  :=  DepartVel ocNew[2]  *  Ml 


Arri veVel ocNew[l]  :=  Arri veVel ocNew[l]  *  M2; 

Arri veVel ocNew[2]  :=  Arri veVel ocNew[2]  *  P2; 

(* 

cal culate  q  *) 

hstar[l]  := 

DeparturePosit01d[2]  *  DepartVel oc01d[3] 
-  DeparturePosit01d[3]  *  DepartVel oc01d[2] 

hstar[2]  := 

DeparturePositOl d[3]  *  DepartVeloc01d[l’J 
-  DeparturePosit01d[l]  *  DepartVeloc01d[3] 

hstar[3]  := 

DeparturePosit01d[l]  *  DepartVel oc01d[2] 
-  DeparturePosit01d[2]  *  DepartVeioc01d[l] 

temp  := 

RlxR2[l]  *  hstar[l] 

+  RlxR2[2]  *  hstar[2] 

+  RlxR2[3]  *  hstar[3]  ; 
if  temp  >=  0.0  then  q  :=  1.0  else  q  :=  -1.0  ; 

end;  (*  GetAl phaBetaGamma  *) 


procedure  UnitVectorh; 
var 

d:  real  ; 


begi  n 

d  :  =  rl  *  r2  *  s i n (AngleBetweenVectors )  ; 

RlxR2[l]  :=  DeparturePosit01d[2]  *  Arri valPosit01d[3] 

-  DeparturePosit01d[3]  *  Arri valPosit01d[2] 

RlxR2[2]  :=  DeparturePosit01d[3]  *  Arri valPosit01d[l] 

-  DeparturePosit01d[l]  *  Arri valPosit01d[3] 


RlxP2[3]  :=  DeparturePosit01d[l]  *  Arri va 1 Posi tOl d [ 2] 
-  DeparturePositQld[2]  *  Arri valPosit01d[l] 


h[l]  :=  RlxR2[l]  /  d  ; 
h [2]  :=  RlxR2[2]  /  d  ; 
h[3]  :=  RlxR.2[3]  /  d  ; 


end;  (*  UnitVectorh  *) 


function  arccos(a: real ):  real; 
var 

temp:  real ; 
begi  n 

if  a  =  U.O  then  a  :=  0.0000000001;  (*  avoid  divide  by  zero  * 

temp  :=  arctan(sq  rt(1.0  -  a  *  a)  /  a); 

if  temp  <  0.0  then  arccos  :=  temp  +  Pi  else  arccos  :=  temp; 
end; 


procedure  AngleRlR2; 
begi  n 

AngleBetween Vectors 


end;  (*  AngleRlR2  *) 


:=  arccos( 

(  DeparturePosit01d[l] 
+  DeparturePosit01d[2] 
+  DeparturePositCl  d[3] 


Arri val Posit01d[l] 
Am  valPosit01d[2] 
Arri val Posit01d[3]  ) 
/  (  rl  *  r2  )  ); 


procedure  UnitVectorcrlr2; 
var 

irl,  ir2,  ic,  c:  array[1..3]  of  real; 
cx  :  rea  1  ; 
temp  :  real ; 
begi  n 

i r  1  [ 1 ]  :=  DeparturePosit01d[l]  /  rl; 
irl [23  :=  DeparturePosit01d[2]  /  rl; 
irl[3]  :=  DeparturePosit01d[3j  /  rl; 


ir2[l]  :=  Arri valPosit01d[l]  /  r2; 
ir2[2]  :=  Am  valPosit01d[2]  /  r2; 
ir2[3]  :=  Arri valPosit01d[3]  /  r2; 


c[l]  :=  Arri valPosit01d[l]  -  0eparturePosit01d[l]; 
c[2]  :=  Arri val Posi t01d[2]  -  DeparturePosit01d[2]; 
c[3]  :=  Arri valPosit01d[3]  -  DeparturePosit01d[3]; 

cx  :=■  sq  rt  (  c[l]  *  c[l]  +  c[2]  *  c[2]  +  c[3]  *  c[3]  ) 
CC  cx; 

i c[ 1 ]  :=  c[l]  /  cx; 
ic[2]  :=  c[2]  /  cx; 
i c[3]  :=  c[3]  /  cx; 

icPirl[l]  :=  ic[l]  +  irl[l]; 
icPirl[2]  :=  ic[2]  +  irl[2]; 
icPirl[3]  :=  i c[3]  +  irl[3]; 


temp  :=  sq  rt(i  cPi  rl  [  1  ]  *  icPiri[l]  +  icPirl[2]  *  icPirl[2] 
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+  i cPi rl[3]  *  icPirl[3]); 


:=  temp; 

IJniti cPi rl[l]  :=  icPirl[l]  /  temp  ; 

Uniti cPirl[2]  :=  icPirl[2]  /  temp  ; 

Uni ti cPi rl[3]  :=  icPirl[3]  /  temp  ; 

icMirl[l]  :=  ic[l]  -  i rl[l] ; 

icMirl[2]  :=  ic[2]  -  irl[2]; 

icMirl[3]  :=  ic[3]  -  irl[3]; 

p  :=  sq rt (i cMi rl[l]*i cMi rl[l]  +  i cMi rl[2]*i cMi rl[2] 

+  i cMi ri[3]*i cMi rl[3]) ; 

Uniti cM i r  1  [ 1 ]  :=  icMirl[l]  /  temp  ; 

Uniti cMi rl[2]  :=  icMirl[2]  /  temp  ; 

IJniti  cMi  rl[3]  :=  icMirl[3]  /  temp  ; 

Ml  :=  temp; 

icPir2[l]  ;=  ic[l]  +  ir2[l]; 

icPir2[2]  :=  i c[2]  +  ir2[2]; 

icPir2[3]  ;=  ic[3]  +  i r 2 [ 3 j ; 

ip  :=  sq  rt (i cPi r2[l]*i cPi r2[l]  +  i cPi r2[2]*i cPi r2[2] 

+  i cPi r2[3]*i cPi r2[3]) ; 

Uniti cPi r2[l]  :=  i cPi r2[ 1 ]  /  temp  ; 

Uniti cPi r2[2]  :=  icPir2[2]  /  temp  ; 

Uniti cPi r2[3]  :=  icPir2[3]  /  temp  ; 


P2  :=  temp; 


i cMi r2[l' 
i  cMi  r2[_2’ 
i cMi r2[3‘ 


ic[l]  -  ir2[l] 
i c[2]  -  i r2[2] 
i c[3]  -  i r2[3] 


ip  :=  sq  rt (i cMi r2[l]*i cMi r2[l]  +  i cMi r2[2]*i cMi r2[2] 

+  i cMi r2[3]*i cMi r2[3]) ; 


Uniti cMi r2[i]  :=  icMir2[l]  /  temp  ; 
Uniti cMi r2[2]  :=  icMir2[2]  /  temp  ; 
Uniti cMi r2[3]  :=  icMir2[3]  /  temp  ; 


M2  :=  temp; 

I;  (*  UnitVectorcrlr2  *) 
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q,|pr.|,  of  M.NOAODS-IM,)-* 


procedure  GetNextSetT rueAnomal ies ; 
var 

counter  :  integer; 
begi  n 

for  counter  :=  1  to  3  do 
begin 

Departu rePos itOld [counter]  :=  OeparturePositi on [TrueStar, counter]; 

Arri valPosit01d[counter]  :=  Arri valPosition[TrueP, counter]; 

DepartVeloc01d[counter]  :=  DepartureVel ocity  [TrueStar, counter]; 

Arri veVel oc01d[counter]  :=  Arri val Velocity [TrueP, counter]; 

end; 

rl  :=  sq  rt(departurepositold[l]  *  departurepositold[l] 
+departurepositold[2]  *  departurepositold[2] 
+departurepositold[3]  *  departurepositold[3]) ; 

r2  :=  sq  rt(arri  valposito1d[i]  *  arri  valpositold[l] 

+arri valpositoid[2]  *  arri valpositold[2] 

♦arri valpositold[3]  *  arri val positold[3]) ; 

end;  (*  proc  GetNextSetT rueAnomal ies  *) 


procedure  FindDeltaV(  AlphaStar,  BetaStar,  GammaStar, 

AlphaP,  BetaP,  GammaP,  C,q,  SG,  rl,  r2.  Theta:  real); 


DeltaVl,  DeltaV2:  real; 
sqrQeltaVl,  sq  rDeltaV2:  real; 
oldk  1  ,oldk 2:  real; 

GravConstStar,  GravConstP:  real; 
GravConstEarth:  real; 

K ,  K 1 ,  K2:  real ; 

Epsilon:  real  ; 

TempEpsilon:  real; 
s :  rea  1 ; 

GravConstSun:  real ; 


function  tan  (a:  real):  real; 
begin 

tan  : =  si n (a  )/cos (a  ) ; 
end ; 


procedure  FindEpsilon; 
var 

accuracy:  real; 

Epsilonl,  Epsilon2:  real; 


procedure  PreEpsilon; 
begin 


Epsilonl  :=  q  /(2.0  *  K  *  (k  1  -  k2))  * 

(  kl  *  AlphaStar  *  (1.0  +  ( r 2  *  cos (theta ) )/C 

+  k2  *  AlphaP  *  (1.0  +  (rl  *  cos (theta ) )/C 

Epsi lon2  :=  SG  /  (2.0  *  K  *  (k  1  -  k2))  * 

(  kl  *  BetaStar  *  (1.0  -  (r2  *  cos (theta ) )/C 

+  k2  *  BetaP  *  (1.0  -  (rl  *  cos (theta ) )/C 

end;  '*  Pre Epsilon  *) 

begin  (*  FindEpsilon  *) 

accuracy  :=  0.00001; 

PreEpsi Ion; 
repeat 

TempEpsilon  :=  ABS(Epsiion)  ; 

Epsilon  :=  arctan(ABS(Epsi1onl  *  sin(TempEpsilon)  + 
write( 

until  (abs (TempEpsilon  -  Epsilon)  <  accuracy); 
end;  (*  FindEpsilon  *) 


begin  (*  main  procedure  FindDeltaV  *) 

TempEpsilon  :=  1.0;  (*  artificial  starting  value  *) 

Epsilon  :=  -0.01;  (*  artificial  starting  value  *) 

GravConstSun  :=  132718E6; 

GravConstEart.h  :=  398603.0; 

GravConstStar  :=  GravConstEarth; 

kl  :=  0.8; 
k 2  :=  0.7; 

s  :=  0.5  *  (rl  +  r2  +  C); 

k  :=  sqrt(  (GravConstSun  *  C)  /  (2.0  *  s  *  (s  -  C)  ) 


-  rl /C  ) 

-  r2/C  )  ); 

+  rl/C  ) 

+  r2/C  )  ); 


Epsilon2) ); 


FindEpsi Ion; 


sqrDeltaVl  :  = 

(K  *  q  /  cos (Epsilon)  -  AlphaStar  ) 

*  (K  *  q  /  cos(Epsilon)  -  AlphaStar  ) 

+  (K  *  SG  *  tan(Epsilon)  -  BetaStar  ) 

*  (K  *  SG  *  tan(Epsilon)  -  BetaStar  ) 

+  (GammaStar  *  GammaStar)  ; 

sqr0eltaV2  :  = 

(K  *  q  /  cos(Epsilon)  -  Alpha?  ) 

*  (K  *  q  /  cos(Epsilon)  -  AlphaP  ) 

+  (K  *  SG  *  tan (Epsilon)  -  BetaP  ) 

*  (K  *  SG  *  tan (Epsilon)  -  BetaP  ) 

+  (GammaP  *  GammaP)  ; 

ol  dk  1  :  =  k  1  ; 

kl  :=  1.0  /  sq  rt(  sqrDeltaVl  +  (2.0  *  Gra  vConstStar  /  6700.0)  ); 
k 2  :=  -1.0  /  sq  re (  sqrDeltaV2  ); 

FindEpsilon; 

until  (ab$(oldkl  -  kl)  <  0.00001); 
sqrDeltaVl  :  = 

(K  *q  /  cos (Epsi  1  on )  -  AlphaStar  ) 

*  (K  * q  /  cos (Epsi Ion)  -  AlphaStar  ) 

/  PI  /  PI 

+  (K  *  SG  *  tan(Epsilon)  -  BetaStar  ) 

*  (K  *  SG  *  tan(Epsilon)  -  BetaStar  ) 

/  Ml  /  Ml 

+  (GammaStar  *  GammaStar)  ; 


sqrDeltaV2  :  = 

(K  *  q  /  cos(Epsilon)  -  AlphaP  ) 

*  (K  * q  /  cos(tpsilon)  -  AlphaP  ) 
/  M2  /  M2 

+  (K  *  SG  *  tan(Epsilon)  -  BetaP  ) 

*  (K  *  SG  *  tan(Epsilon)  -  BetaP  ) 

/  P2  /  P2 

+  (GammaP  *  GammaP)  ; 


TrueStar*10:4,TrueP*10:4,9q  rt(sq  r0eltaVl):12:2,9q  rt(sq  rOeltaV2):l2:2) 


writeln( 'END' ); 

end;  (*  proc  FindDeltaVl  ■*) 


procedure  DeltaV360; 
var 

counter:  integer; 
begin 

TrueStar  :=  11  ; 
write( 'S6  =  ’ ) ; 
readln(sg); 

for  counter  :=  0  to  360  do 
begin 

TrueP  :=  counter; 

GetNextSetT  rueAnoma lies; 

UnitVectorcrlr2; 

AngleRlR2; 

UnitVectorh; 

GetAl phaBetaGamma ; 

FindDe1taV( 

Depart Vel ocNew[l], Depart Vel ocNew[2],DepartVel ocNew[3], 
Arri veVelocNew[l],Arri veVelocNew[2],Arri veVel ocNew[3], 
CC  *  AU^},SG,rl  *  AU,r2  *  AU .AngleBetweenVectors ) ; 

end;  (*  for  counter  *) 

end;  (*  proc  DeltaV350  *) 


begin  (*  main  *) 
InitialDisplay ; 

GecDataFi le; 

SelectObject ; 
ConvertToRadians; 
PreComputeT  ranscendentals ; 
Departure360; 
InitialDisplay : 

GetDataFi le ; 

SelectObject; 
ConvertToRadians ; 
PreComputeT  ranscendentals ; 
Arri val 360; 

Cl rScr; 

DeltaV360; 


writeln; 

writeln; 

readln;  (*  return  to  DOS  *) 


end.  (*  main  *) 
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Appendix  C. 


PASCAL  Program  to  Generate 
the  Va 1 i dati on  Ephemerides 

program  Hel i oCentri c0rbit2050; 


(*  This  program  calculates  ephemerides  for  solar  system  objects.  It  *) 
(*  uses  the  following  data  files:  PLANET. ELE,  COMET. ELE,  AST-0001. ELE  *) 
(*  AAA82-84.ELE  and  any  user-defined  data  file  with  the  same  format.  *) 
(*  Either  a  single  date  or  a  period  ephemeris  can  be  output.  *) 


const 

MaxRecSize  =  25; 

Obi  iq  Eel  ip  =  23.439; 
rAdjustl985  =  -5479.5; 
JD00Janl985  =  2446065.5; 


type 

Name  =  string  [10]; 
ElementSet  =  record 

AsteroidName:  Name; 

Tp:  real; 
i  :  real; 
w  :  real ; 
N  :  real; 
a  :  real ; 
e  :  real ; 
Mo:  real; 

end; 


(*  object  name  *) 

(*  date  of  perihelion  passage /epoch  *) 
(*  inclination  to  ecliptic  *) 

(*  arguement  of  perihelion  *) 

(*  longitude  of  ascending  node  *) 

(*  semi-major  axis,  in  AU  *) 

(*  ec.centri city  *) 

(*  mean  ancmoly  at  To, 

Mo  =  0  at  perihelion  *) 


var 

AsteroidFi 1 e :  text; 

FileName:  string[14]; 

ElementSetRec:  array  [1. .MaxRecSize]  of  ElementSet; 
Ifor:  integer; 


MaxNoOfAsteroids:  i 

nteger; 

Selected:  integer; 

(*  record  selected  to  work  with  *) 

A1,B1,Y1,  A2,B2,Y2: 

real 

(* 

auxiliary  quantities  *) 

Px.Py.Pz,  Qx,Qy,Qz: 

real 

f  * 

\ 

Gaussian  Constants  *) 

M  : 

real 

(* 

Mean  anomoly  *) 

no: 

real 

\ 

Mean  daily  mot ion, degrees  per  day  * 

Eccentri cAnomaly : 

real 

i* 

tc  solve  Kepler's  Equation  *) 

X ref sun, Y ref sun, Z ref sun: 

real 

(* 

i  n  AU  * ) 

Rrefsun: 

real 

(* 

'•adius  from  sun,  in  AU  *) 

Xsun  ,Ysun,Zsun: 

real 

l* 

i  n  AU  * ) 

Xref geo,Yref geo.Zrefgeo: 

real 

(* 

i  n  AU  * ) 

R ref geo: 

rea  1 

i  * 

\ 

sun  to  earth  distance  in  AU  *) 

RA,Dec: 

real 

(* 

right  ascention,  declination  *) 

Jul ianOate: 

real 

AsteroidName: 

name; 

Tp: 

real ; 

(* 

date  of  perihelion  passage/epoch  *) 

rea  1 ; 
real ; 
real ; 
real ; 
rea  1 ; 


Mo:  real; 


(*  inclination  to  ecliptic  *) 

(*  arguement  of  perihelion  *) 

(*  longitude  of  ascending  node  *) 

(*  semi-major  axis,  in  AU  *) 

(*  eccentricity  *) 

(*  mean  anomoly  at  To, 

Mo  =  0  at  perihelion  *) 


function  int(x:real):  integer;  (*  int  corresponds  to  BASIC  ' i nt 1  *) 
begin 

if  (x  >=  0)  or  (x  =  trunc(x))  then  int  :=  trunc(x) 
else  int  :=  trunc(x)  -  1 
end;  (*  int  *) 


procedure  GetDate(var  Year, Mon, Day  , Hour: real ) ; 


reply  =  string[3]; 

month  =  (Jan, Feb, Mar, Apr, May  .Jun.Jul , Aug, Sep, Oct, Nov, Dec, 

NotYetAssigned); 

const  Blank  =  -1;  (*  for  dumny  variables  *) 

MinYear  =  1900;  (*  limits  of  Julian  Date  formula  *) 
MaxYear  =  2099;  (*  limits  of  Julian  Date  formula  *) 

var 

WhatMonth:  reply; 

MaxDays:  integer; 

StartMonth,  EndMonth:month; 

StartDay  .EndDay ,StartHour,EndHour: integer; 

GetYear,  GetDay :  integer; 

GetMonth:  month; 

GetHour:  real; 

Days.PhaseDays:  real; 


begi  n 


GetYear  :=  Blank  ; 

while  not  (GetYear  in  [MinYear. .MaxYear])  do 
begin 

write( 'What  year,  '.MinYear,'  ...  '.MaxYear,'  ?  '); 

readln(GetYear); 
end;  (*  whi le  *) 

Year  :=  GetYear; 

while  GetMonth  =  NotYetAssigned  do 
begi  n 

write('What  month,  Jan. .Dec  ?  '); 

readln (WhatMonth ) ; 

if  WhatMonth  =  'Jan'  then  GetMonth  :=  Jan; 

if  WhatMonth  =  'Feb'  then  Getmonth  :=  Feb; 

if  WhatMonth  =  'Mar'  then  Getmonth  :=  Mar; 

if  WhatMonth  =  'Apr'  then  Getmonth  :=  Apr; 


if  WhatMonth  =  'May'  then  GetMonth 
if  Whatmonth  =  'Jun'  then  GetMonth 
if  WhatMonth  =  * Ju 1 1  then  GetMonth 
if  WhatMonth  =  'Aug'  then  GetMonth 
if  WhatMonth  =  ‘Sep1  then  GetMonth 
if  WhatMonth  =  'Oct'  then  GetMonth 
if  WhatMonth  =  'Nov'  then  GetMonth 
if  WhatMonth  =  ‘Dec*  then  Getmonth 
end;  (*  while  *) 
case  GetMonth  of 

Jan, Mar, May  ,Jul , Aug, Oct .Dec:  MaxDays  :  = 

Sep, Apr, Jun, Nov:  MaxDays  :=  30; 

Feb:  MaxDays  :=  29 
end;  (*  case  *) 


51; 


May ; 
Jun; 
Jul ; 
Aug; 
Sep; 
Oct; 
Nov; 
Dec 


case  GetMonth  of 

Jan:  Mon 

=  1.0 

Feb:  Mon 

=  2.0 

Mar:  Mon 

=  3.0 

Apr:  Mon 

=  4.0 

May :  Mon 

=  5.0 

Jun:  Mon 

=  6.0 

Jul:  Mon 

=  7.0 

Aug:  Mon 

=  8.0 

Sep:  Mon 

=  9.0 

Oct:  Mon 

=  10.0 

Nov:  Mon 

=  11.0 

Dec:  Mon 

=  12.0 

end;  (*  case  *) 

GetDay  : =  Blank ; 

while  not  (GetDay  in  [1. .MaxDays])  do 
begi  n 

write('What  day,  1.. '  .MaxDays, '  ?  '); 

readl n (GetDay  ) 
end;  (*  wni le  *) 

Day  : =  GetDay ; 

GetHour  :  =  B lank  ; 

while  not  ((GetHour  >=  0  )  and  (GetHour  <=  24.0)  )  do 
begi  n 

write(‘What  nour,  0.0  ...  24.0  ?  '); 

readln(GetHour) 
end;  (*  whi le  *) 

Hour  :=  GetHour; 
end;  (*  GetDate  *) 


procedure  CalendarDate(JulianDate:real ;  var  Day  .Month, Year: integer) 

(*  worfcs  for  dates  1100-2000  *) 

(*  appears  slightly  inaccurate  before  1583  *) 

label  1,2  ; 
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JulianDate  :=  JulianDate  -2400000.0  ; 

NDate  :=  JulianOate  -  15018.0  +0.0001  ;  {*  .0001  fcr  roundoff  error  *) 
Yearl  :=  NDate/365.25  ; 

Year  :=  1900  +  int(Yearl)  ; 


if  (Year/4-int(Year/4)=0  )  and  (Year/100-i nt (y ear/100)<>0  )  then 

LeapYear  :=  1; 

NDay  :=  365.25  *(Yearl  -  int(Yearl)  }  ; 

if  NDay  -  int(N0ay)>  0.5  then  NDay  :=  NDay  +  1.0  ; 

Day  :  =  int(NDay)  ; 

if  LeapYear  =  0  then  Day  :=  Day  -  1  ; 

Day  :=  Day  +  int((Year-2000)/100)  ; 

if  Year  <  1583  then  Day  :=  Day  -  (10+int((Year  -  1500)  /  100))  ; 
if  Day  -  32  <  0  then  begin 

Montn  :=  1  ; 

Day  :=  Day  -  0  ; 

goto  1  ; 
end; 

if  Day  -  60  <  0  then  begin 

Month  :=  2  ; 

Day  Day  -  31  ; 
goto  1  ; 
end; 

if  Day  -  91  <  0  then  begin 

Month  :=  3  ; 

Day  :=  Day  -  59  ; 
goto  1  ; 
end; 

if  Day  -  121  <  0  then  begin 

Month  :=  4  ; 

Day  :=  Day  -  90  ; 
goto  1  ; 
end; 

if  Day  -  152  <  0  then  begin 

Month  :=  5  ; 

Day  :=  Day  -  120  ; 
goto  1  ; 
end; 

if  Day  -  182  <  0  then  begin 

Month  :=  6  ; 

Day  :=  Day  -  151  ; 
goto  1  ; 
end; 

if  Day  -  213  <  0  then  begin 
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Month  :=  7  ; 

Day  :=  Day  -  181  ; 
goto  1  ; 
end; 

if  Day  -  244  <  0  then  begin 

Month  :=  8  ; 

Day  :=  Day  -  212  ; 
goto  1  ; 
end; 

if  Day  -  274  <  0  then  begin 

Month  :=  9  ; 

Day  :=  Day  -  243  ; 
goto  1  ; 
end; 

if  Day  -  305  <  0  then  begin 

Month  :=  10  ; 

Day  :=  Day  -  273  ; 
goto  1  ; 
end; 

if  Day  -  335  <  0  then  begin 

Month  :=  11  ; 

Day  :  =  Day  -  304  ; 
goto  1  ; 
end; 

if  Day  -  366  <  0  then  begin 

Month  :=  12  ; 

Day  :=  Day  -  334  ; 
end; 

1 : i f  (LeapYear  =  1)  and  (Month  >  2)  then  Day  :=  Day  -  -  ; 

if  Day  >=  1  then  goto  2  ;  ^ 

if  (Day  <  1)  and  (LeapYear  =  1)  and  (Month  =  o)  then  begin  ^ 

Day  •  ~  ^  > 
Month  :=  2  ; 
goto  2  ; 
end; 

if  (Day  <  1)  and  (LeapYear  =  0)  and  (Month  =  3)  then  begin  ^  _ 

uay  *  -  , 

Month  :=  2  ; 
goto  2  ; 
end; 

if  (day  <  1)  and  (Month  =  1)  then  begin 

Day  :=  31  ; 

Month  :=  12  ; 

Year  :=  Year  -  1  ; 
goto  2  ; 
end; 

if  (day  <  1)  and  ((Month  =  2)  or  (Month  =  4)  or  (Month  =  6) 
or  (Month  =  9)  or  (Month  =11))  then  begin 

Day  :  =  31  ; 

Month  :=  Month  -  1; 
goto  2  ; 
end; 
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if  Day  <  1  then  begin  - 

Day  :=  30  ; 

Month  :=  Month  -  1  ; 
end; 

2:  if  Day  =  365  then  begin 

Day  :=  31  ; 

Month  :=  12  ; 
end; 

if  Day  =  366  then  begin 

Day  :=  1  ; 

Month  :=  1  ; 

Year  :=  Year  +  1  ; 
end; 

end;  (*  CalendarOate  *) 


procedure  GetJul i anDate  (Year, Month, Day ,Uni vTime: real ; 

var  JulianDate: real ) ; 

(*  input:  Universal  Time  and  Date  from  28  Feb  1900  to  31  Dec  2099  *) 

(*  output:  Julian  Date  2  400  000  .  0000  *) 

begin  (*  GetJul ianDate  *) 

JulianDate  :=  367  *  Year  -  trunc(7/4  *  (Year  +  trunc((Month  +  9)  /  12))) 

+  trunc(275  *  Month  /  9)  +  Day  +  1721013.5 
+  UnivTime  /  24  ; 


end;  (*  GetJulianDate  *) 


procedure  Sol veKeplensEq  n  (M:real;  var  Eccentri cAnomaly  : real ) ; 
const 

accuracy  =  0.0000000001; 
var 

Mtemp:  real; 

Mrad:  real; 

Iteration:  real; 

diff:  real;  (*  temporary  difference  *) 

Eccentri cAnomalyRadians:  real; 

begin 

Mtemp  :=  -1.0;  (*  artificial  starting  value  *) 

Mrad  :>  M*Pi/180.0; 

Eccentri cAnomalyRadians  :=  Mrad  -  0.1  ;  (*  rough  starting  point  *) 


wnile  Mtemp  <  Mrad  do 


begi  n 

Eccentri  cAnomalyRadians  :=  Eccentri  cAnomalyRadians  +  0.1; 

Mtemp  :=  Eccentri cAnomaly Radians  -  e*sin(Eccentri cAnomalyRadians ); 
end; 

while  Mtemp  >  Mrad  do 
begi  n 

Eccentri cAnomalyRadians  :=  Eccentri cAnomalyRadians  -  0.1; 

Mtemp  :=  Eccentri cAnomalyRadians  -  e*sin(Eccentri cAnomalyRadians ) ; 
end; 

iteration  :=  0.1; 

diff  :=  1.0;  (*  artificial  start  *) 
while  diff  >  accuracy  ao 
begi  n 

iteration  :=  iteration/2.0; 
if  Mtemp  <  Mrad 
then 

Eccentri cAnomalyRadians  :=  Eccentri cAnomalyRadians  +  iteration 
else 

Eccentri cAnomalyRadians  :=  Eccentri cAnomalyRadians  -  iteration 

Mtemp  :=  Eccentri cAnomalyRadians  -  e*sin(Eccentri cAnomalyRadians ) ; 
diff  :=  abs(Mrad  -  Mtemp); 
end;  (*  while  *) 

Eccentri cAnomaly  :=  Eccentri cAnomalyRadians  *  180.0/Pi; 
end;  (*  procedure  Sol veKeplersEq n  *) 


procedure  SunXYZcoord(Jul ianDate:  real;  var  Xcoord:  real; 

var  Y coord:  real; 
var  Z coord:  real); 
var 

L :  rea  1 ; 
g:  real; 
gamma :  rea  1  ; 

E :  rea  1 ; 

R :  rea  1 ; 
n:  real  ; 

begi  n 

n  :=  nAdjustl985  +  JulianOate  -  JD00Janl985; 
g  :=  357.528  +  0.9856003  *  n; 

while  g  <  0.0 

do  g  :=  g  +  360.0  ; 
while  g  >  360.0 
do  g  g  -  360.0  ; 

L  :=  280.460  +  0.9856474  *  n  ; 


while  L  <  0.0 

do  L  :=  L  +■  360.0  ; 
while  L  >  360.0 
do  L  ;=  L  -  360.0  ; 

gamma  :=  L  +  1.915  *  sin(g*Pi /180.0)  +  0.020  *  sin(2.0  *  g  *  Pi/180.0); 
E  :=  23.439  -  0.0000004  *  n  ; 

R  :=  1.00014  -  0.01671  *  cos(g  *  Pi/180.0)  -  0.00014  *  cos(2.0  *  g 

*  Pi/180.0); 

Xcoord  :=  R  *  cos (gamma  *  Pi/180.0)  ; 

Ycoord  :=  R  *  cos(E  *  Pi/180.0)  *  sin(gamma  *  Pi/180.0)  ; 

Zcoord  :=  R  *  sin(E  *  Pi/180.0)  *  sin(gamma  *  Pi/180.0)  ; 


end  ;  (*  SunXYZcoord  *) 


procedure  GetDataFile;  (*  Reads  data  into  global  variables  *) 

(*  ***  Same  as  procedure  GetDataFile  **********  *) 

(*  ***  in  program  Hulk owerLauBenderDeltaVee  ***  *) 


□rccedure  SelectOoject ;  (*  Global  variables  *) 

/*  ***  came  as  procedure  SelectObject  *********  *) 
(*  ***  in  program  Hulk owerLauBenderDel taVee  ***  *) 


procedure  Cal culateObjPosition  (OateToCal culate:  real); 
begi  n 

SunXYZcoord (DateT  oCal culate,Xsun, Ysun,Zsun ) ; 

E ccentri cAnoma ly  :=  0.0; 

A1  :=  sin(N*Pi /180.0)*sin(w,tPi  /180.0);  (*  auxi  1  iary  q uantities  *) 

B1  :=  cos (N*Pi /180.0)*sin(w*Pi /180.0) ; 

Y1  :=  si n (i *Pi /180.0)*si n (w*Pi /180.0) ; 

A2  sin (N*Pi /180.0)*cos (w*Pi /180.0) ; 

B2  :=  cos (N*Pi /180.0)*cos (w*Pi /ISO. 0 ) ; 
v2  :=  sin (i *Pi /180.0)*cos (w*Pi /180.0) ; 

(*  GAUSSIAN  CONSTANTS  *) 

Px  :=  B2  -  Al*cos (i *Pi /180.0) ; 

Py  :=  (A2  +  Bl*cos (i *Pi /180.0) )*cos (Obi iq Eel i p*Pi /180.0)  - 

Yl*sin(0bliqEclip*Pi  /180.0); 

Pz  :=  (A2  +  Bl*cos (i *Pi /18Q.0) )*sin{0bl iq  Eel i p*Pi  /ISO . 0 )  + 

Yl*cos(0bliqEclip*Pi/180.0); 


k\' 


:=  -B1  -  A2*cos (i *Pi /180.0) ; 

:=  (-A1  +  82*cos (i *Pi /180.0) )*cos (ObliqEcl ip*Pi /180.0)  - 

Y2*sin(Obl  iq  Eel  ip*Pi  /1 80.0) ; 

:=  (-A1  +  B2*cos(i*Pi/180.0))*sin(0bliqEclip*Pi/180.0)  + 

Y2*cos (Obi  iq  Eel i p*Pi /180.0 ) ; 


no  :=  0.9856076686/sq  rt  (a*a*a  ) ;  (*  mean  daily  motion  in  degrees  *) 

M  :=  Mo  +  no*(DateToCa1  dilate  -  Tp);  (*  mean  anomaly,  DateToCalculate  is 

time  of  ephemeris  *) 

Sol  veKeplersEq  n  (M,  Eccentri  cAnomaly  ) ; 

Xrefsun  :=  a*Px*(cos (Eccentri cAnomaly *Pi /180.0)-e )  + 

a*sq  rt  ( l-e*e )*Qx*si n (Eccentri cAnomaly  *Pi /180.0) ; 

Yrefsun  :=  a*Py *(cos (Eccentri cAnona ly *Pi /180. 0 )-e )+ 

a*sq  rt  ( l-e*e  )*Qy  *s  in  (Eccentri  cAnomaly  *Pi /180.0) ; 

Zrefsun  :=  a*Pz *(cos (Eccentri cAnoma’y *Pi /180.0)-e )  + 

a*sq  rt (l-e*e )*Qz *s in (Eccentri cAnomaly *Pi /180.0) ; 

Rrefsun  :=  sq  rt  (Xrefsun*Xrefsun  +  Yrefsun*Yrefsun  +  Zrefsun*Zrefsun ) ; 

(*  in  AU  *) 


Xrefgeo  :=  Xrefsun  +  Xsun  ; 

Yrefgeo  :=  Yrefsun  +  Ysun  ; 

Zrefgeo  :=  Zrefsun  +  Zsun  ; 

Rrefgeo  :=  sq  rt. (Xrefgeo*Xref geo  +  Yrefgeo*Yref geo  *  Zref geo*Zref geo) ; 

(*  in  AU  *) 

RA  :=  arctan(Yref geo/Xref geo)*180.0/Pi ;  (*  RA  in  degrees  *) 

if  Xrefgeo<0  then  RA  :=  RA  +  180.0 

else  if  Yrefgeo<0  then  RA  :=  RA  +  360.0  ;  (*  correct  q uadrant  *) 

RA  RA/ 1 5.0;  (*  RA  in  hours  *) 

DEC  :=  arctan( (Zrefgeo/Rref geo)/sq  rt(1.0  -  Zref geo/Rrefgeo*Zref geo 

/Rrefgeo))  *  180.0/Pi;  (*  in  degrees  *) 

end;  (*  CalculateObjPos ition  *) 


procedure  InitialDisplay ; 

(*  ***  Same  as  procedure  InitialDisplay  *******  *) 
(*  ***  in  program  Hulk owerLauBenderDel taVee  ***  *) 


procedure  Peri odEphemeris ; 
var 

Year, Month, Day .Hour:  real; 
Ephemerislnterval :  real; 
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DateToCal culate:  real; 

StartJD.EndJD:  real; 

Choice:  char; 

PrinterOrScreen:  char; 

RAhour,  RAmin,  DECdeg,  DECmin,  DECsec:  integer; 
RAsec:  real; 


procedure  Print(DateToPrint:  real); 
var 

Day .Month, year:  integer; 
begin  (*  Print  *) 

CalendarDate(DateToPrint .Day .Month .Year ) ; 

if  Pri nterOrScreen  in  ['S', 's'] 
then 
begin 

write  (Day : 2, Month: 3, Year: 5, ' 
writeln(RAhour:3,RAmin:3,RAsec:5:l,DECdeg:6, 

DECmi n: 3, DECsec: 3  ); 

end  (*  then  *) 
else 
begi  n 

write(LST,Day : 2, Month: 3, Year: 5, ' 

wri tel n (LST .RAhour : 3 .RArni n: 3 .RAsec: 5: 1 .DECdeg: 6, 

DECmin:3,DECsec:3  ); 

end;  (*  else  *) 
end;  (*  Print  *) 


procedure  Cal cEphemeris ; 

procedure  TruncateRADEC; 
begin 

RAhour  :=  trunc(RA); 

RAmin  :=  trunc(f rac(RA)*60) ; 

RAsec  :=  frac(f rac(RA)*60)*60; 

DECdeg  :=  trunc(DEC); 

DECmin  :=  abs (trunc(frac(DEC )*60) ) ; 

DECsec  :=  abs (trunc(f rac(f rac(DEC )*60)*60) ) ; 
end;(*  TruncateRADEC  *) 

begin  (*  Cal cEphemeris  *) 

DateToCalculate  :=  StartJD; 

while  DateToCalculate  <  EndJD  do 
begin 

Cal culateObjPositi on (DateToCal culate); 
TruncateRADEC; 


Print (Da teToCa leu  late); 

DateToCal culate  :=  DateToCal culate  +  Ephemerislnterval; 
end;  (*  whi le  *) 

DateToCal culate  :=  EndJD; 

Cal cul ateOb jPositi on (Da teToCa lcul ate ) ; 

TruncateRADEC; 

Print (DateT  oCal culate ) ; 

end;  (*  Cal cEphemeris  *) 


begin  (*  PeriodEphemeris  *) 
writeln; 

writeln( 'Starting  Date'); 

Get Da te( Year .Month , Day  .Hour ) ; 

GetJul ianDate( Year, Month .Day  .Hour.Jul ianDate) ; 
StartJD  :=  JulianDate; 
repeat 
writeln; 

writeln('For  a  SINGLE  position,  type  <S> ' ) ; 
writeln('For  MULTIPLE  positions,  type  <M>‘); 
readln(Choi  ce) 

until  Choice  in  [ 'S 1 , 's 1 , ‘M* , 'm ']; 

if  Choice  in  ['S', 's']  then 
begi  n 

EndJD  :=  StartJD; 

Epnemerislnterval  :=  1.0;  (*  Dumny  *) 


begi  n 
wri  tel  n ; 

writeln(  'End  Date' ); 

GetDate(Year, Month, Day .Hour); 

GetJul ianDate (Yea r, Month, Day  .Hour ,Jul ianDate) ; 

EndJD  :=  JulianDate; 
writeln; 

write('What  is  interval  for  ephemeris  in  decimal  days?  '); 
readln (Ephemerislnterval ); 
end;  (*  else  *) 

writeln; 

repeat 

wri tel n( 'Output  to  <S>  Screen  or  <P>  Printer  ??'); 
readl n (Pri nterOrScreen ) ; 
until  PrinterOrScreen  in  ['S ' , 's ' , * P * , 'p ' ]; 

if  PrinterOrScreen  in  ['S', 's']  then 
begin 

writeln; 


writeln( 'Ephemeris  for  >>>  .AsteroidName ,  <<<  ); 

writeln('DD  MM  YEAR  Right  Ascen  Declination'); 


end 


else 


begi  n 

writeln(LST ) ; 

writeln(LST, ‘Ephemeris  for 
writeln(LST, ' DD  MM  YEAR 
end; 


>>>  ‘  ,AsteroidName, 1  «<*) 
Right  Ascen  Declination') 


Cal cEphemeris; 
end;  (*  Peri odEphemeris  *) 


begin  f*  main  *) 
InitialDisplay ; 

Get.DataFi  le; 

SelectObject ; 

Peri odEphemeris ; 

writeln; 

writeln; 

readln;  (*  return  to  00S  *) 
end.  (*  main  *) 


Appendix  D. 


PASCAL  Program  to  Generate  True  Anomal ies  verses  Date 
program  Cal cAnomal ies ; 


(*  Calcs  True  Anomalies  every  10  days  from  Oh  -  05  Jan  1985  to  Oh  *) 
(*  -  2  Jan  2020  for  Earth  and  Target.  Stores  values  in  Vearth[0. .1278]*) 
(*  and  Vtarget[0..1278].  Saves  them  in  user-designated  file  *) 
(*  Note:  element  set  must  be  entered  directly  in  code.  Program  does  *) 
(*  not  read  .ELE  files.  *) 


const 

usun  =  1.327 18E20 ;  (*  Grav  const  of  sun  in  m3/s2  *) 

All  =  1 . 49600E II;  (*  Astronaumi  cal  Unit  in  m  *) 

var 

Mearth  :  real;  (*  Mean  anomaly  earth  at  epoch,  in  degrees  *) 

Tearth  :  real;  (*  Time  of  epocn  earth  *) 

eearth  :  real;  (*  Eccentri  city  earth  *) 

aearth  :  real;  (*  Semi-major  axis  a  of  earth  in  m  *) 

Vearth  :  array [0. .2000]  of  real;  (*  True  anomaly  of  earth  *) 

Mtarget  :  real;  (*  Mean  anomaly  target  at  epoch,  in  degrees  *) 

Ttarget  :  real;  (*  Time  of  epoch  target  *) 

etarget  :  real;  (*  Eccentri  city  target  *) 

atarget  :  real;  (*  Semi-major  axis  a  of  target  in  m  *) 

Vtarget  :  array [0. .2000]  of  real;  (*  True  anomaly  of  target  *) 

counter  :  integer; 

JD  :  real;  (*  Julian  date  *) 

Name  :  string[8]; 


function  tan (x: real  ): real  ; 
begin 

tan  :=  sin(x)/cos(x); 
end; 


procedure  KeplerEq  uati on(M,e: real ,  var  X:real); 

(*  Adapted  from  : 

Danby,  J.M.A  and  T.M.  Burkardt.  "The  Solution  of  Kepler's  Equation,  I," 
Celestial  Mechanics,  31:  95-107  (May  1983).  *) 

var 


ES: 

real ; 

F: 

real ; 

EC: 

real ; 

FP: 

real ; 

OX: 

real ; 

(*  Initial  Guess  *) 


writer  ~  ) ; 

X  :=  M; 

ES  :=  e  *  sin(X) ; 

F  :=  X  -  ES  -  M; 

repeat 

EC  :=  e  *  cos(X); 

FP  :=  1.0  -  EC; 

DX  :  =  -F/FP; 

DX  :=  — F / (FP  +  DX*ES/2.0); 

OX  :=  -F/(FP  +  DX*ES/2.0  +  DX*DX*EC/6.0) ; 

DX  :=  -  F  /  (FP  +  DX  *  ES/2.0  +  OX  *  DX  *  EC/6.0 

-  DX  *  DX  *  DX  *  ES/24.0) ; 

X  :=  X  +  DX; 

ES  :=  e  *  sin(X) ; 

F  :=  X  -  tS  -  M; 
until  ABS(f)  <=  0.0000000001; 

end ; 


procedure  Bacon(M,T,e,a:  real;  z:  integer); 
var 


counter  : 

integer; 

tp:  real; 

(*  Time  since  perihelion  passage 

MX :  real  ; 

(*  Calculated  mean  anomaly  *) 

X  :  rea  1 ; 

(*  Eccentric  anomaly  *) 

V  :  real  ; 

(*  True  anomaly  *) 

Too:  real  ; 

(*  Time  of  perihelion  passage  * 

n  :  real;  (*  Mean  motion  *) 

JD  :  real;  (*  Current  Julian  date  for  calculations  *) 


begi  n 


n  :=  sqrt  (  usun  /  (a*a*a)  )  *  3600.0  *  24.0  ; 

tp  :=  M  /  n  ; 

Tpo  :=  T  -  tp; 

JD  :=  2446066.5;  (*  Oh  -  01  Jan  1985  *) 

JD  :=  2446066.5  +  4.0;  (*  Start  on  mult  of  10  (2446070.5) 

Oh  -  05  Jan  1985  *) 


counter  :=  0; 


while  Tpo  >  JD  do 
begin 

Tpo  :=  Tpo  -  2.0  *  Pi  /  n; 
end; 


while  JD  <  2458850.5  do 


(*  Oh  -  02  Jan  2020  *) 


MX  :=  n  *  (JO  -  Tpo); 


while  MX  >=  (2.0  *  Pi )  do 
begi  n 

MX  :=  MX  -  2.0  *  Pi; 

Tpo  :*  Tpo  +  2.0  *  Pi  /  n; 
end; 

KeplerEquation(MX,e,X); 

V  :=  2.0  *  arctan(tan(0.5  *  X)  * 

sq  rt(  (1.0  +  e  )/(1.0  -  e)  )  ); 

if  V  <  0.0  then  V  :=  V  +  2.0  *  Pi  ; 

V  :=  V*180.0/Pi  ; 

if  z  =  1  then  Vearth[ccunter]  :=  V  else  Vtarget[counter]  :=  V; 
JD  :=  JD  +  10.0; 

counter  :=  counter  +  1; 

end;  (*  while  JD  <  244...  *) 
end;  (*  procedure  Bacon  *) 


procedure  SaveAnomaly Fi le; 
const 

MaxNoAnomal ies  =  1278  ; 
type 

AnomRec  =  record 

Ve  :  rea  1 ; 

Vt  :  real ; 
end; 
var 

AnomFile  :  file  of  AnomRec; 

Anom  :  AnomRec; 
counter  :  integer; 
begi  n 
writeln; 

writeln( 'Input  prefix  for  anomaly  data  file 
readln(Name); 

assi gn(AnomF i le,Name+‘ .DTA* ) ; 

Rewri  t.e(AnomF  i  le ) ; 
with  Anom  do 
begi  n 

for  counter  :=  0  to  MaxNoAnomal ies  do 
begin 

Ve  :=  Vearth[counter] ; 

Vt  :=  Vtarget[counterj; 
write (AnomFile, Anom); 
end; 

end;  (*  with  Anom  *) 
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cl ose (AnomFi 1 e ) ; 

end;  (*  Procedure  SaveAnomalyFi le  *) 


procedure  LoadAnomalyFile; 
const 

MaxNoAnomalies  =  1278  ; 
type 

AnomRec  =  record 

Ve  :  real ; 

Vt  :  real ; 
end; 
var 

AnomFile  :  file  of  AnomRec; 

Anom  ;  AnomRec; 
counter  :  integer; 
begi  n 

assign(AnomFile, 'ANOM-l.DTA' ) ; 

Reset (AnomFi le ) ; 
with  Anom  do 
begi  n 

for  counter  :=  0  to  MaxNoAnomalies  do 
begin 

read (AnomFi le,Anom) ; 
Vearth[counter]  :=  Ve; 
Vtarget[counter]  :=  Vt; 
end; 

end;  (*  with  Anom  *) 
close(AnomFile); 

end;  (*  Procedure  LoadAnomalyFile  *) 


begin 

Mearth  :=  208.89818  *  Pi/180.0; 

Tearth  :  =  2446280.5;  (*  All  1985  Astronomical  Almanac  data  *) 

eearth  :=  0.0166912; 
aearth  :=  1.496G328E11; 

bacon(Mearth,  Tearth,  eearth, aearth, 1); 

(*  **  Asteroid  1982  XB  **  Hulk  ewer,  Lau,  Bender  data  *) 

Mtarget  :=  266.0197  *  Pi /ISO. 0; 

Ttarget  ;=  2446000.5; 
etarget  :=  0.4468762; 
atarget  :=  1.837667  *  AU; 

bacon (Mtarget,  Ttarget,  etarget, atarget, 0); 

JD  :=  2446070.5; 
writeln; 


(*  Remove  brlaces  here  for  simultaneous  listing  *) 

(*  for  counter  :=  0  to  1278  do  *) 

(*  2446070.5  +  10.0  *  1278  =  2458850.5  *) 

(*  begin 

writeln(JD  +  10.0*counter: 12: 1 ,Vearth[counter]:8: 1 , 

Vtarget [counter]: 8: 1 ) ; 

end;  *) 

SaveANoma’yFile; 


9} U  }  Ji.iJI  1»  I  .1.  1  ■■U»l  ■  f 


Appendix  E. 

PASCAL  Program  to  List  True  Anomal ies  verses  Date 
program  ListAnomal  ies; 

(*  List  True  Anom  every  10  days  1985-2020.  05  Jan  1985  -  02  Jan  2020  *) 
(*  Lodes  for  file  AN0M-1.DTA.  Lists  True  Anomalies  for  both  departure  ) 
(*  and  arrival  planets  in  separate  tabular  form.  ) 

var 

Vearth,  Vtarget,  V:  array [0.. 1278]  of  real; 


procedure  LoadAnomalyFile; 
const 

MaxNoAnomal ies  =  1278  ; 
type 

AnomRec  =  record 

Ve  :  real ; 

Vt  :  real ; 
end; 
var 

AnomFile  :  file  of  AnomRec; 

Anom  :  AnomRec; 
counter  :  integer; 

begin 

assi gn(AnomFi le, 'AN0M-1.DTA ' ) ; 

Reset (AnomFile); 
with  Anom  do 
begi  n 

for  counter  :=  0  to  MaxNoAnomal ies  do 
begi  n 

read(AnomFi le,Anom) ; 
Vearth[counter]  :=  Ve; 
VtargetCcounter]  :=  Vt; 
end; 

end;  (*  with  Anom  *) 
cl ose(AnomFi le) ; 

end;  (*  Procedure  LoadAnomalyFile  *) 


procedure  ListAnom; 
const 

JD  =  2446070.5; 
var 

counter,  column  :  integer; 
begi  n 

for  counter  :=  0  to  127  do  (*  2446070.5  +  100.0  *  127  =  2458770.5  *) 
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if  (  counter  =  0)  or  (counter  =  43)  or  (counter  =  86)  then 

begi  n 

readln; 

wri tel n (LST ) ;  writeln(LST) ; 
writeln(LST);  writeln(LST);  writeln(LST); 

writeln(LST,  _  .... 

•  Julian  Date  70  80  90  00  10  20  30  40  50  60') 

writeln(LST, 

I  _  _ _ _ _ _ _ _ _ _  -  -  ) 

end; 

write (LST, '  '  ,JD  +  100.0*counter: 12: 1 , 1  '); 

for  column  :=  0  to  9  do 
begi  n 

write(LST,V[counter+co1umn]:4:0) ; 
end; 

wri  tel  n  (LST ) ; 

end; 

end;  (*  procedure  ListAnom  *) 


begi  n 

LoaoAnomalyFi  le; 

V  :=  dearth; 
ListAnom; 

V  :=  Vtarget; 
ListAnom; 


Appendix  F. 

PASCAL  Program  to  Find  Launch  Dates 
program  FindLaunchDate; 

(*  Finds  a  match  of  2  true  anomalies  and  time-of-f 1 i ght  if  a  match  *) 
(*  exists  from  January  1985  to  January  2020.  Resolution  is  10  degrees  *) 
(*  of  true  anomaly  and  10  day  intervals  of  dates.  Each  search  is  for  *) 

(*  true  anomaly  of  departure  plus/minus  10  degrees,  true  anomaly  of  *) 

(*  arrival  plus/minus  10  degrees,  and  time-of-f 1 i ght  plus  or  minus  10  *) 
(*  days.  A  total  of  27  combinations  are  searched  (3x3x3)  for  each  set  *) 
(*  of  true  anomalies  and  T0F  entered.  Uses  data  file  AN0M-1.DTA  or  *) 

(*  user  supplied  data  file.  This  data  file  contains  1278  entries  of  *) 

(*  true  anomalies  every  10  days  from  Jan  1985  to  Jan  2020  for  both  *) 
(*  departure  and  arrival  planets.  *) 

const 

MaxNoAnomal ies  =  1278  ; 

MaxArraySize  =  1279  ; 

var 

LaunchFile  :  text; 

LaunchFileName  :  string[12]; 

DeltaV  :  real; 

M  :  string[9]; 

D,Y  :  integer; 

TAEarth  :  array [1. .MaxArraySize]  of  integer; 

TATarget  :  array [1. .MaxArraySize]  of  integer; 

TAE,TAT,  ITAEarth,  ITATarget  :  integer; 

T0F,  T0F1 ,  T0F2  :  integer; 

Agn  :  cnar  ; 

Again  :  boolean; 

procedure  LoadAnomalyFi le; 
type 

AnomRec  =  record 

Ve  :  real  ; 

Vt  :  real; 
end; 

var 

AnomFile  :  file  of  AnomRec; 

Anom  :  AnomRec; 
counter  :  integer; 

AnomFileName  :  string[12]; 

begi  n 

writeln;  writeln( 'Input  prefix  of  Anomaly  file  '); 

readl n (AnomF i leName ) ; 

ass i gn (AnomF i 1 e , AnomF i 1 eName+ 1 .DTA 1 ) ; 

Reset (AnomF i le ) ; 


begi  n 

for  counter  :=  1  to  MaxNoAnomal ies  do 
begi  n 

rea d { AnomF i le.Anom) ; 

TAEarth[counter]  :=  round(Ve/10.0) ; 

TATarget[counter]  :=  round{Vt/10.0) ; 
end; 

end;  (*  with  Anom  *) 
close(AnomFile); 

end;  (*  Procedure  LoadAnomalyFi le  *) 
procedure  JulianToDate(JD:real  ); 

(*  Convert  Julian  dates  to  calendar  dates.  *) 

(*  Works  from  1  Jan  1900  to  near  end  of  21st  century.  *) 

(*  Writes  date  in  form:  e.g.  12  September  1984  *) 

(*  Julian  Date  procedures  adapted  from  public  domain  software:  author  *) 

(*  unknown  *) 

var  date  :  real ; 

Jul ian, Day  .Month, Year:  integer; 


procedure  JtoD(Julian:  integ er;var  Day  .Month, Year:  integer); 

(*  Convert  from  a  Julian  date  to  a  calenda-  date  *) 

(*  Note  that  much  care  is  taken  to  avoid  problems  with  inaccurate  bit 
representations  inherent  in  the  binary  fractions  of  the  real  numbers 
used  as  temporary  variables.  Thus  the  seemingly  unnecessary  use  of 
small  fractional  offsets  and  int()  functions  *) 

var  Temp:  real ; 
begi  n 

Temp  :=  i nt (32767. 5+Jul ian ) ;  (*  Convert  16  bit  quantity  to  real  no.  *) 
if  Temp <58. 5 
then 

begin  (*  The  first  two  months  of  the  twentieth  century  are 

handled  as  a  special  case  of  the  general  algorithm 
used  which  handles  all  of  the  rest  *) 

Year  :=  1900; 
if  Temp <30. 5 
then 
begin 

Month  :=  1; 

Day  :=  round(Temp+1.0) 
end 
else 
begi  n 

Month  :=  2; 

Day  :=  round(Temp-30.0) 


end 
el  se 
begi  n 

Teinp  :=  i  nt  (4.0*(Temp-59.0)+3.5) ; 

Year  :=  trunc(Temp/1461. 0+0. 00034223) ; 

(*  0.00034223  is  about  one  half  of  the  reciprocal  of  1461.0  *) 

Day  :=  succ(round(Temp-Year*1461.0)  div  4); 

Month  :=  (5*Day-3)  div  153; 

Day  :=  succ( (5*Day-3)  mod  153  div  5); 

Year  :=  Year+1900; 
if  Month<10 
then 

Month  :=  Month+3 
else 
begin 

Month  :=  Month-9; 

Year  :=  succ(Year) 
end 
end 
end ; 


procedure  WriteDate(Julian:  integer); 

(*  Write  the  date  out  to  the  console  in  long  form  , 
e.g.  "10  September  1984"  *) 

const 

Months:  array [1. . 12]  of  string[9]  = 

( ‘January  'February  V  March', 'April  ', 'May  ',  'June ' , 

'July ' , 'August ' , 'September' , 'October' , 

'November' , 'December' ); 
var  Day .Month, Year  :  integer; 

begi  n 

JtoD(Julian, Day .Month, Year);  (*  Convert  into  date  form  *) 

write(Day,'  ' ,Months[Month], '  '.Year); 

D  :  =  Day  ; 

M  :=  Months [Month]; 

Y  :=  Year; 

end ; 


begin  (*  JulianToDate  *) 

JD  :=  JD  -  2447787.99999; 
Julian  :=  round(jd); 

Wri teDate(Jul ian ) ; 


end;  (*  JulianToDate  *) 


procedure  InputTAandTOF; 
begin 

TATarget[1279]  :=  1000;  (*  Dumny  for  search  routine  *) 

writeln; 

write( ‘Input  True  Anomaly  of  Departure  Planet  in  degrees  '); 
readln(ITAEarth); 

ITAEarth  :=  round{ITAEarth/10) ; 
writeln; 

write( 'Input  True  Anomaly  of  Arrival  Planet  in  degrees  '); 
readl  n  (ITATarger. ) ; 

ITATarget  :=  round(ITATarget/10) ; 
writeln; 

write('Input  Delta  V  in  km/sec  '); 
readln (DeltaV) ; 
wr i tain ; 

write('Input  Time-of-Fl i ght  in  days  ’); 

readl n (TOF ) ; 

writeln; 


TOF  :=  round(T0F/10) ;  (*  Time-of-Fl i ght  in  10*days  *) 

T0F1  :=  TOF  +  1; 

T0F2  :=  TOF  -  1; 
end;  (*  InputTAandTOF  *) 

procedure  DoCombos; 
var 

counter  :  integer; 
procedure  Check  Combos; 

procedure  WriteDetai ls(increment:  integer); 
begi  n 

JulianToDate(2446070.5  +  10.0*counter); 
wri tel n (TAE*10:4,TAT*10:4 , 

I0,r(T0F+increment ) :  5,  TOF :  6 ,  DeltaV:4:l); 

writeln (LaunchFile,D: 2, '  ',M,‘  ' , Y : 4 , ‘  *  ,TAE*10:4,TAT*10:4, 

10*(T0F+increment ): 5,  T0F:6,  De 1 1 a V : 4 : 1 ) ; 
end;  (*  WriteDetai Is  *) 


procedure  ChedcTATarget; 
begi  n 

if  TATarget[counterxTOF] 

1 f  TATarget[counter+TOF+l] 
if  TATarget[counter+TOF-l] 
end;  (*  CheckTATarget  *) 


ITATarget  then  WriteDetai Is (0) ; 
ITATarget  then  WriteDetai Is (1 ) ; 
ITATarget  then  WriteDetai Is (-1 ) ; 


begin  (*  CheckCombos  *) 
for  TAE  :=  ITAEarth-1  to  ITAEarth+l  do 


•-v-*Cy 

NV-VV 

% *-  JC 


egi  n 

if  TAE  =  -1  then  TAE  :=  35; 

if  TAE  =  37  then  TAE  :=  10; 

for  TAT  :=  ITATarget-1  to  ITATarget+1  do 
begin 

if  TAT  =  -1  then  TAT  :=  35; 

if  TAT  =  37  then  TAT  :=  10; 

if  TAEarth[counter]  =  ITAEarth  then  ChedcTATarget 
end;  (*  for  TAT  *) 
end;  (*  forTAE  *) 
end;  (*  CheckCombos  *) 

begin  (*  DoCombos  *) 
counter  :=  1; 

while  counter  <  1279-TOF  do 
begi  n 

CheckCombos ; 
counter  :=  counter  +  1; 
end;  (*  while  counter  x) 
end;  (*  DoCombos  *) 


procedure  InitLaunchFi le; 

begin 

writeln; 

writeln( 'Input  prefix  of  Launch  file  '); 
readln(LaunchFileName-); 

assign(LaunchFi le,LaunchFi leName+'.LAU 1 ) ; 
Rewri te(LaunchFi le) ; 
end;  (*  Procedure  InitLaunchFile  *) 


begin  (*  Main  FindLaunchDate  *) 
LoadAnomalyFi  le; 

InitLaunchFi le; 
again  :=  true; 
whi  le  Again  =  true  do 
begi  n 

InputTAandTOF ; 

DoCombos ; 
writeln; 

writeln( 'Again  ?  (Y/N) 

readln(Agn); 

if  upcase(agn)  =  ' N '  then 
end; 

close(LaunchFile); 
wri teln ( 'END  END'); 
read! n; 

end.  (*  Main  FindLaunchDate  *) 


?  (Y/N)  '); 

=  ' N '  then  Again  :=  false; 
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224 

207 

190 

174 

158 

142 

127 

112 

97 

300 

286 

266 

247 

228 

210 

191 

174 

156 

140 

124 

108 

93 

290 

310 

283 

260 

237 

216 

196 

176 

158 

141 

124 

109 

96 

280 

347 

310 

279 

252 

228 

205 

185 

166 

148 

133 

121 

112 

270 

385 

340 

304 

273 

246 

222 

201 

183 

167 

156 

143 

118 

260 

424 

375 

336 

302 

274 

250 

229 

214 

202 

179 

152 

132 

250 

465 

416 

374 

341 

313 

290 

273 

262 

229 

198 

176 

163 

240 

508 

460 

420 

388 

363 

344 

333 

296 

260 

235 

220 

215 

230 

550 

507 

471 

444 

423 

410 

385 

342 

313 

297 

292 

627 

220 

585 

548 

519 

497 

484 

479 

447 

413 

395 

616 

688 

771 

210 

603 

575 

554 

540 

534 

538 

532 

510 

667 

727 

796 

866 

200 

596 

577 

564 

559 

562 

575 

596 

685 

733 

788 

842 

893 

190 

562 

550 

545 

548 

558 

578 

663 

701 

743 

786 

827 

860 

180 

503 

498 

500 

508 

526 

608 

638 

670 

705 

737 

765 

787 

170 

428 

430 

437 

452 

533 

556 

582 

609 

636 

659 

679 

692 

160 

353 

359 

372 

451 

470 

491 

514 

537 

558 

575 

589 

596 

150 

286 

297 

423 

419 

411 

429 

449 

469 

486 

499 

508 

512 

i40 

233 

391 

387 

379 

368 

378 

396 

413 

427 

437 

442 

442 

130 

359 

356 

348 

337 

323 

341 

357 

372 

382 

339 

391 

388 

120 

327 

320 

309 

296 

300 

317 

331 

342 

349 

352 

351 

348 

no 

295 

284 

272 

273 

289 

303 

313 

320 

323 

322 

319 

314 

100 

263 

251 

255 

271 

283 

291 

297 

298 

297 

294 

290 

283 

90 

232 

246 

260 

270 

275 

278 

278 

276 

272 

267 

261 

253 

80 

242 

254 

260 

262 

263 

261 

259 

255 

250 

243 

235 

226 

70 

250 

253 

253 

252 

249 

246 

242 

236 

230 

222 

214 

204 

60 

251 

249 

246 

243 

238 

233 

228 

221 

214 

206 

196 

137 

50 

250 

245 

240 

225 

230 

224 

217 

209 

201 

192 

183 

169 

40 

250 

244 

227 

231 

224 

217 

208 

200 

191 

179 

169 

158 

30 

253 

245 

236 

228 

220 

211 

202 

192 

182 

171 

160 

149 

20 

258 

247 

238 

228 

217 

207 

197 

186 

175 

164 

153 

141 

10 

264 

252 

240 

228 

216 

204 

193 

181 

169 

157 

146 

134 

0 

273 

257 

243 

229 

215 

201 

189 

176 

164 

152 

140 

128 

Arr'^’Dep 

240 

250 

260 

270 

280 

290 

300 

310 

320 

330 

340 

350 

360 

116 

105 

94 

83 

73 

64 

55 

48 

44 

45 

43 

184 

350 

no 

98 

87 

75 

65 

55 

47 

42 

42 

40 

188 

183 

340 

103 

91 

79 

68 

57 

48 

41 

39 

36 

193 

189 

259 

330 

96 

84 

71 

60 

49 

41 

37 

32 

201 

211 

295 

370 

320 

89 

76 

54 

52 

43 

38 

31 

211 

241 

346 

432 

486 

310 

83 

70 

58 

48 

42 

32 

221 

277 

406 

510 

572 

594 

300 

80 

68 

60 

48 

38 

231 

309 

461 

588 

667 

697 

692 

290 

85 

76 

59 

49 

241 

334 

497 

647 

753 

804 

811 

784 

280 

94 

76 

65 

251 

357 

521 

684 

814 

892 

926 

910 

865 

270 

100 

88 

286 

389 

543 

707 

849 

949 

1009 

1024 

996 

940 

260 

120 

116 

434 

573 

729 

868 

979 

1058 

1101 

1103 

1066 

1005 

250 

158 

491 

616 

757 

389 

995 

1078 

1136 

1164 

1158 

1116 

1053 

240 

558 

669 

794 

912 

1005 

1081 

1133 

1182 

1201 

1188 

1145 

1086 

230 

725 

831 

930 

1013 

1075 

1121 

1164 

1195 

1204 

1190 

1150 

1094 

220 

858 

942 

1009 

1060 

1093 

1125 

1153 

1174 

1179 

1163 

1128 

1077 

210 

931 

986 

1026 

1053 

1070 

1087 

1108 

1123 

1123 

1108 

1077 

1032 

200 

937 

967 

989 

1001 

1006 

1016 

1030 

1041 

1040 

1024 

996 

959 

190 

886 

904 

911 

913 

913 

918 

929 

935 

933 

917 

893 

863 

180 

801 

808 

808 

804 

799 

802 

809 

815 

810 

796 

774 

747 

170 

699 

700 

696 

688 

681 

681 

688 

690 

685 

671 

654 

629 

160 

599 

595 

587 

577 

568 

568 

574 

575 

568 

559 

538 

516 

150 

510 

504 

493 

480 

469 

470 

475 

475 

473 

456 

437 

418 

140 

438 

429 

417 

402 

388 

391 

396 

401 

387 

371 

354 

337 

130 

383 

373 

360 

343 

323 

331 

345 

334 

319 

303 

287 

272 

120 

341 

331 

320 

303 

274 

303 

293 

278 

263 

249 

235 

221 

no 

308 

300 

292 

284 

234 

248 

241 

230 

218 

206 

194 

182 

100 

276 

267 

254 

225 

203 

204 

200 

193 

183 

173 

162 

150 

90 

243 

232 

208 

192 

178 

175 

171 

164 

156 

146 

136 

125 

80 

215 

204 

183 

169 

158 

154 

148 

142 

133 

124 

115 

105 

70 

194 

176 

163 

152 

142 

136 

130 

123 

115 

107 

98 

89 

60 

171 

160 

148 

138 

129 

122 

115 

108 

100 

92 

84 

76 

50 

158 

147 

136 

126 

117 

no 

102 

95 

87 

80 

72 

65 

40 

147 

1  37 

126 

116 

107 

99 

91 

84 

76 

69 

62 

56 

30 

138 

128 

117 

107 

98 

89 

81 

74 

66 

60 

54 

51 

20 

130 

120 

109 

99 

89 

80 

72 

65 

58 

52 

49 

50 

10 

123 

112 

101 

91 

81 

72 

63 

56 

50 

47 

48 

47 

0 

116 

105 

94 

33 

73 

64 

55 

48 

44 

45 

43 

184 

Appendix  J 


True  Anomaly  verses  Date  from  the  Year 
1985  to  2020  for  the  Asteroid  1985TB 


Julian  Date 

70 

80 

90 

00 

10 

20 

30 

40 

50 

60 

2446070.5 

217 

218 

220 

221 

222 

224 

225 

227 

228 

230 

2446170.5 

232 

234 

236 

238 

240 

242 

244 

247 

249 

252 

2446270.5 

255 

258 

262 

265 

269 

274 

279 

284 

289 

296 

2446370.5 

303 

310 

319 

328 

337 

347 

358 

8 

19 

28 

2446470.5 

38 

46 

54 

61 

68 

74 

79 

84 

89 

93 

2446570.5 

97 

100 

104 

107 

110 

112 

115 

117 

119 

122 

2446670.5 

124 

126 

127 

129 

131 

132 

134 

136 

137 

138 

2446770.5 

140 

141 

142 

144 

145 

146 

147 

148 

149 

150 

2446870.5 

152 

153 

154 

155 

156 

157 

157 

158 

159 

160 

2446970.5 

161 

162 

163 

164 

165 

165 

166 

167 

168 

169 

2447070.5 

170 

170 

171 

172 

173 

174 

174 

175 

176 

177 

2447170.5 

177 

178 

179 

180 

181 

181 

182 

183 

184 

184 

2447270.5 

185 

186 

187 

188 

188 

189 

190 

191 

192 

192 

2447370.5 

193 

194 

195 

196 

157 

197 

198 

199 

200 

201 

2447470.5 

202 

203 

204 

205 

206 

207 

208 

209 

210 

211 

2447570.5 

212 

213 

214 

216 

217 

218 

219 

221 

222 

224 

2447670.5 

225 

227 

228 

230 

232 

233 

235 

237 

239 

241 

2447770.5 

244 

246 

249 

252 

254 

258 

261 

265 

269 

273 

2447870.5 

278 

283 

288 

295 

301 

309 

317 

326 

336 

346 

2447970.5 

356 

6 

17 

27 

36 

45 

53 

60 

67 

73 

2448070.5 

78 

84 

88 

92 

96 

100 

103 

106 

109 

112 

2448170.5 

114 

117 

119 

121 

123 

125 

127 

129 

131 

132 

2448270.5 

134 

135 

137 

138 

140 

141 

142 

143 

145 

146 

2448370.5 

147 

148 

149 

150 

151 

152 

153 

154 

155 

156 

2448470.5 

157 

158 

159 

160 

161 

162 

163 

164 

164 

165 

2448570.5 

166 

167 

168 

169 

169 

170 

171 

172 

173 

173 

2448670.5 

174 

175 

176 

177 

177 

178 

179 

180 

180 

181 

2448770.5 

182 

183 

184 

184 

185 

186 

187 

187 

188 

189 

2448870.5 

190 

191 

191 

192 

193 

194 

195 

196 

196 

197 

2448970.5 

198 

199 

200 

201 

202 

203 

204 

205 

206 

207 

2449070.5 

208 

209 

210 

211 

212 

213 

214 

215 

217 

218 

2449170.5 

219 

221 

222 

223 

225 

226 

228 

230 

231 

233 

2449270.5 

235 

237 

239 

241 

243 

246 

248 

251 

254 

257 

2449370.5 

260 

264 

268 

272 

277 

282 

287 

293 

300 

308 

2449470.5 

316 

324 

334 

344 

354 

5 

15 

25 

34 

43 

2449570.5 

51 

59 

66 

72 

78 

83 

87 

92 

96 

99 

2^49670.5 

103 

106 

109 

111 

114 

116 

119 

121 

123 

125 

2449770.5 

127 

129 

130 

132 

133 

135 

136 

138 

139 

141 

2449870.5 

142 

143 

144 

146 

147 

148 

149 

150 

151 

152 

2449970.5 

153 

154 

155 

156 

157 

158 

159 

160 

161 

162 

2450070.5 

163 

163 

164 

165 

166 

167 

168 

168 

169 

170 

2450170.5 

171 

172 

172 

173 

174 

175 

176 

176 

177 

178 

2450270.5 

179 

179 

180 

181 

182 

183 

183 

184 

185 

186 

Julian  Date  70  80  90  00  10  20  30  40  50  60 


2450370.5 

2450470.5 

2450570.5 

2450670.5 

2450770.5 

2450870.5 

2450970.5 

2451070.5 

2451170.5 

2451270.5 

2451370.5 

2451470.5 

2451570.5 

2451670.5 

2451770.5 

2451870.5 

2451970.5 

2452070.5 

2452170.5 

2452270.5 

2452370.5 

2452470.5 

2452570.5 

2452670.5 

2452770.5 

2452870.5 

2452970.5 

2453070.5 

2453170.5 

2453270.5 

2453370.5 

2453470.5 

2453570.5 

2453670.5 

2453770.5 

2453870.5 

2453970.5 

2454070.5 

2454170.5 

2454270.5 

2454370.5 

2454470.5 

2454570.5 


186  187 

195  195 

204  204 
214  215 
228  229 
248  251 
286  292 

13  23 
87  91 
118  120 
136  138 
149  150 
159  160 
167  168 
175  176 
183  184 
191  192 
200  201 
209  210 
221  223 
238  240 
266  271 
330  340 
63  70 
108  110 
130  131 
144  145 
155  156 
164  165 
172  173 
180  131 
188  189 

196  197 

205  206 
216  217 
230  232 
252  255 
296  303 

29  39 
93  97 
122  124 
139  140 


188  189 
196  197 
205  206 
216  218 
231  233 
253  256 
299  306 

33  42 

95  99 
123  125 
139  140 
151  152 
161  162 
169  170 
177  178 
185  186 
193  194 
201  202 
212  213 
224  226 
243  245 
275  280 
350  1 

76  81 
113  116 
133  134 
146  147 
157  158 
166  167 
174  175 
182  182 

189  190 
198  199 
207  208 
219  220 
234  236 
259  262 
311  320 

47  55 
101  104 
126  128 
141  143 


190  190 

198  199 
207  209 
219  220 
235  237 
260  263 
314  323 

50  58 
102  105 
126  128 
142  143 
153  154 
162  163 
171  172 
179  179 
186  187 
194  195 
203  204 
214  215 
227  229 
247  250 
285  291 

11  21 

86  90 
118  120 
136  137 
149  150 
159  160 
167  168 
175  176 
1S3  184 

191  192 

199  200 
209  210 
221  223 
238  240 
266  270 
329  338 

62  69 
107  110 
129  131 
144  145 


151  152  153  154  155  156 


191  192 
200  201 
210  211 
222  223 
239  241 
267  271 
332  342 
65  71 

108  111 
130  132 
144  145 
155  156 
164  165 
172  173 
180  181 
188  189 
196  197 
205  206 
216  217 
231  232 
253  256 
298  305 
31  40 

94  98 
122  124 
139  140 
151  152 
160  161 
169  170 
177  178 
135  185 
193  193 
201  202 
211  212 
224  226 
242  244 
274  279 
349  359 
75  80 

113  115 
133  134 
146  147 
157  158 


193  194 

202  203 
212  213 
225  226 
243  245 
276  281 
352  3 

77  82 
113  116 
133  135 

147  148 
157  158 
166  167 
174  175 
182  182 
190  190 
198  199 
207  208 
219  220 
234  236 
259  263 
313  321 

49  56 
101  105 
126  128 
141  143 
153  154 
162  163 
171  171 
178  179 
186  187 

194  195 

203  204 
214  215 
227  229 
247  250 
284  290 

9  20 
85  89 
117  120 
136  137 

148  150 
159  159 


Julian  Date 

70 

80 

90 

00 

10 

20 

30 

40 

50 

60 

2454670.5 

160 

161 

162 

163 

164 

165 

166 

166 

167 

168 

2454770.5 

169 

170 

170 

171 

172 

173 

174 

174 

175 

176 

2454870.5 

177 

178 

178 

179 

180 

181 

181 

182 

183 

184 

2454970.5 

185 

185 

186 

187 

188 

188 

189 

190 

191 

192 

2455070.5 

192 

193 

194 

195 

196 

197 

198 

198 

199 

200 

2455170.5 

201 

202 

203 

204 

205 

206 

207 

208 

209 

210 

2455270.5 

211 

212 

213 

215 

216 

217 

218 

220 

221 

222 

2455370.5 

224 

225 

227 

228 

230 

232 

234 

235 

237 

240 

2455470.5 

242 

244 

246 

249 

252 

255 

258 

261 

265 

269 

2455570.5 

273 

278 

283 

289 

295 

302 

310 

318 

327 

337 

2455670.5 

347 

357 

8 

18 

28 

37 

46 

54 

61 

68 

2455770.5 

74 

79 

84 

39 

93 

97 

100 

103 

107 

109 

2455870.5 

112 

115 

117 

119 

121 

123 

125 

127 

129 

131 

2455970.5 

132 

134 

135 

137 

138 

140 

141 

142 

144 

145 

2456070.5 

146 

147 

148 

149 

150 

151 

152 

154 

155 

155 

2456170.5 

156 

157 

158 

159 

160 

161 

162 

163 

164 

165 

2456270.5 

165 

166 

167 

168 

169 

170 

170 

171 

172 

173 

2456370.5 

173 

174 

175 

176 

177 

177 

178 

179 

180 

180 

2456470.5 

181 

182 

183 

134 

184 

185 

186 

187 

188 

188 

2456570.5 

189 

190 

191 

192 

192 

193 

194 

195 

196 

197 

2456670.5 

197 

198 

199 

200 

201 

202 

203 

204 

205 

206 

2456770.5 

207 

208 

209 

210 

211 

212 

213 

214 

216 

217 

2456870.5 

218 

219 

221 

222 

223 

225 

226 

228 

230 

231 

2456970.5 

233 

235 

237 

239 

241 

244 

246 

249 

251 

254 

2457070.5 

257 

261 

264 

268 

273 

277 

282 

288 

294 

301 

2457170.5 

308 

317 

325 

335 

345 

355 

6 

16 

26 

35 

2457270.5 

44 

52 

60 

66 

73 

78 

83 

88 

92 

96 

2457370.5 

100 

103 

106 

109 

112 

114 

117 

119 

121 

123 

2457470.5 

125 

127 

129 

130 

132 

134 

135 

137 

138 

139 

2457570.5 

141 

142 

143 

145 

146 

147 

148 

149 

150 

151 

2457670.5 

152 

153 

154 

155 

156 

157 

158 

159 

160 

161 

2457770.5 

162 

163 

164 

164 

165 

166 

167 

168 

169 

169 

2457870.5 

170 

171 

172 

173 

173 

174 

175 

176 

176 

177 

2457970.5 

178 

179 

180 

180 

181 

182 

183 

183 

184 

185 

2458070.5 

186 

187 

187 

188 

189 

190 

191 

191 

192 

193 

2458170.5 

194 

195 

196 

196 

197 

198 

199 

200 

201 

202 

2458270.5 

203 

204 

205 

206 

207 

208 

209 

210 

211 

212 

2458370.5 

213 

214 

215 

217 

218 

219 

220 

222 

223 

225 

2^58470.5 

226 

228 

229 

231 

233 

235 

237 

239 

241 

243 

2458570.5 

246 

248 

251 

254 

257 

260 

264 

268 

272 

276 

2458670.5 

281 

287 

293 

300 

307 

315 

324 

333 

343 

353 

2458770.5 

4 

14 

24 

34 

43 

51 

58 

65 

iC 


Appendix  K 


T rue  Anomaly  verses  Date  from  the  Year 
1985  to  2020  for  the  Earth 


Julian  Date 

70 

80 

90 

00 

10 

20 

30 

40 

50 

60 

2446070.5 

2 

12 

22 

33 

43 

53 

63 

73 

83 

93 

2446170.5 

102 

112 

122 

131 

141 

151 

160 

170 

179 

189 

2446270.5 

198 

208 

218 

227 

237 

247 

256 

266 

276 

286 

2446370.5 

296 

306 

316 

326 

336 

346 

357 

7 

17 

27 

2446470.5 

37 

47 

57 

67 

77 

87 

97 

107 

117 

126 

2446570.5 

136 

146 

155 

165 

174 

184 

193 

203 

213 

222 

2446670.5 

232 

241 

251 

261 

271 

281 

291 

301 

311 

321 

2446770.5 

331 

341 

351 

1 

12 

22 

32 

42 

52 

62 

24^6870.5 

72 

32 

92 

102 

112 

121 

131 

141 

150 

160 

2446970.5 

169 

179 

138 

198 

207 

217 

227 

236 

246 

256 

2447070.5 

266 

275 

285 

295 

305 

315 

326 

336 
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15 

26 

36 
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46 

56 

66 
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96 
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10 
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40 

51 

61 

71 

2449070.5 

31 

90 
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2449170.5 
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196 

206 

216 

225 

235 

244 
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2449270.5 

274 

284 

294 

304 

314 

324 

334 

344 

354 

5 

2449370.5 

15 

25 

35 

45 

55 

65 

75 

35 

95 

105 

2449470.5 

115 

124 

134 

144 

153 

163 

172 

182 

191 

201 

2449570.5 
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249 

259 

269 

279 
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2449670.5 

309 

319 

329 

339 
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359 

9 

20 

30 

40 
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60 

70 

80 

90 
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263 
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323 
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2450070.5 

344 
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4 

14 

24 

35 

45 

55 

65 

75 
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85 

95 
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162 

172 

2450270.5 
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268 
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99 
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224 
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253 
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283 
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303 
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313 

323 

333 

343 

353 

4 

14 

24 

34 

44 

2450870.5 

54 

64 

74 

84 

94 
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114 
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133 

143 
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162 

171 

181 
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200 

209 

219 

229 

238 

2451070.5 

248 
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287 

297 

307 

318 

328 

338 
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348 
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8 

19 

29 

39 

49 

59 

69 

79 

2451270.5 
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99 

108 

118 

128 

137 

147 

157 

166 
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2451370.5 
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195 

204 

214 

224 

233 

243 

253 

262 

272 

2451470.5 

282 

292 

302 

312 

322 

332 

343 

353 

3 

13 

2451570.5 

23 

33 

44 

54 

64 

74 

84 

93 
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132 

142 
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209 
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228 

238 

247 
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156 

166 

175 

185 

194 

204 

213 

223 

233 

242 

2452170.5 
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93 

103 

113 

122 

132 

142 

151 

161 

170 

180 

2452470.5 
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Comparison  of  Computed  Ephemerides  with  Almanac  Data 
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Appendix  M. 

Orbital  Elements  of  Selected  Solar  System  Objects 


Argue. 

Long. 

Semi  - 

of 

Ascend. 

Major 

Eccen- 

Mean 

Name 

Epoch 

Incl in. 

Perigee 

Node 

Axis 

tri  ci  ty 

Anomaly 

. 

(+2440000. 1 

1  (deg) 

(^9) 

.(dej) 

(AU) 

(deg) 

1985  JA 

6200.5000 

36.73616 

288.7836 

232.012 

1.64345 

0.32013 

288.784 

1985  PA 

6300.5000 

55.92239 

311.4754 

147.371 

1.42205 

0.30137 

263.249 

1985  TB 

6432.5546 

27.02800 

66.8280 

23.390 

2.61136 

0.57289 

000.000 

1982  DB 

5647.7650 

1.42009 

157.8281 

314.082 

1.48932 

0.36017 

0.0000 

1982  XB 

6000.5000 

3.87314 

16.6825 

74.5784 

1.33767 

0.44688 

266.0197 

1943  Ant. 

6400.5000 

8.70253 

338.1039 

245. 785 

1.43019 

0.25591 

128.210 

1982  B6 

6400.5000 

20.94324 

253.6499 

129.2863 

1.40682 

0.35475 

182.2862 

1982  DV 

6000.5000 

5.92683 

349.1949 

218.2229 

2.03324 

0.45667 

315.4827 

1982  FT 

6000.5000 

20.38330 

234.5135 

348.3263 

1.77421 

0.28380 

16.9277 

1982  RA 

6400.5000 

32.97512 

53.2491 

339.4558 

1.57466 

0.28372 

194.8538 

1982  RB 

6000.5000 

24.99579 

158.5766 

158.4470 

2.10193 

0.39487 

263.0551 

1982  TA 

6000.5000 

12.11626 

118.5972 

10.0439 

2.30298 

0.77104 

187.4833 

1982  YA 

6000.5000 

34.57320 

143.6389 

269.1622 

3.70673 

0.69725 

98.2681 

1983  LB 

6000.5000 

25.39914 

220.1517 

80.9363 

2.29142 

0.47869 

128.4299 

1983  LC 

6000.5000 

1.51906 

184.6915 

159.0668 

2.63152 

0.70933 

101.2684 

1983  RB 

6000.5000 

19.42719 

114.8082 

168.8844 

2.22334 

0.50700 

141.7119 

1983  RD 

6000.5000 

9.51734 

192.9485 

173.4031 

2.09011 

0.48667 

129.0814 

1983  SA 

6000.5000 

30.77923 

316.6039 

350.0285 

4.23072 

0.71457 

97.2714 

1983  TB 

6400.5000 

22.02894 

321.6679 

265.0462 

1.27132 

0.89017 

205.4340 

1983  TF 

5620.5000 

7.83661 

121.0523 

10.4301 

1.34276 

0.38708 

283.8069 

1983  VA 

6400.5000 

16.23778 

11.6840 

76.8703 

2.61068 

0.69170 

167.1080 

1984  KB 

6000.5000 

4.63662 

334.8782 

170.5624 

2.22103 

0.76228 

61.6257 

1984  KD 

6400.5000 

13.61579 

203.5824 

81.8423 

2.19758 

0.54087 

155.5356 

1984  QA 

6400.5000 

9.91826 

54.8267 

152.0450 

0.98963 

0.46838 

181.1960 

1  Mercury 

6280.5000 

7.00578 

29.0864 

48.3492 

0.38710 

0.20564 

230.6966 

2  Venus 

6280.5000 

3.39475 

54.8432 

76.7188 

0.72332 

0.00681 

256.0015 

3  Earth 

6280.5000 

0.00185 

107.9279 

354.9000 

1.00002 

0.01669 

208.8982 

4  Mars 

6280.5000 

1.85078 

286.3834 

49.6025 

1.52372 

0.09331 

140.6931 

5  Jupiter 

6280.5000 

1.30465 

275.1670 

100.4667 

5.20266 

0.04808 

301.3302 

6  Saturn 

6280.5000 

2.48456 

339.0884 

113.7135 

9.55775 

0.05121 

141.1483 

7  Uranus 

6280.5000 

0.77457 

102.2384 

74.0564 

19.27850 

0.04662 

75.3400 

8  Neptune 

6280.5000 

1.76881 

230.14C8 

131.8112 

30.25704 

0.00752 

271.57664 

9  Pluto 

6280.5000 

17.13148 

114.1751 

110.4183 

39.60047 

0.25128 

353.66898 

Sources : 

:  1985JA 

(3:1) 

1985PA 

(3:2) 

Planets 

( 20: E3 ) 

1985TB 

(3:3) 

Others 

(14:4) 
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